options(warn = -1)
options(repr.plot.width=15, repr.plot.height=15)

1 Предварительный анализ данных

library(readxl)
col_I_sn<-read_excel("stat2021_sm\\us_col.std\\I_shortname.xls")
## New names:
## • `` -> `...1`
col_I<-read_excel("stat2021_sm\\us_col.std\\I.xls")
## New names:
## • `` -> `...1`
head(col_I_sn)|>rmarkdown::paged_table()
summary(col_I_sn)
##      ...1               PPIND            FICE          STATE          
##  Length:176         Min.   :1.000   Min.   : 1009   Length:176        
##  Class :character   1st Qu.:1.000   1st Qu.: 1738   Class :character  
##  Mode  :character   Median :1.000   Median : 2550   Mode  :character  
##                     Mean   :1.301   Mean   : 2901                     
##                     3rd Qu.:2.000   3rd Qu.: 3417                     
##                     Max.   :2.000   Max.   :10366                     
##                                                                       
##      TYPE              AVRMATH         AVRVERB         AVRCOMB    
##  Length:176         Min.   :390.0   Min.   :391.0   Min.   : 810  
##  Class :character   1st Qu.:520.0   1st Qu.:454.5   1st Qu.: 980  
##  Mode  :character   Median :544.0   Median :479.5   Median :1017  
##                     Mean   :563.5   Mean   :494.7   Mean   :1058  
##                     3rd Qu.:603.0   3rd Qu.:518.5   3rd Qu.:1126  
##                     Max.   :750.0   Max.   :665.0   Max.   :1410  
##                     NA's   :68      NA's   :68      NA's   :67    
##     AVR_ACT          MATH_1          MATH_3          VERB_1     
##  Min.   :19.00   Min.   :350.0   Min.   :480.0   Min.   :320.0  
##  1st Qu.:22.00   1st Qu.:460.0   1st Qu.:590.0   1st Qu.:400.0  
##  Median :23.00   Median :500.0   Median :627.5   Median :430.0  
##  Mean   :23.69   Mean   :515.4   Mean   :635.0   Mean   :448.3  
##  3rd Qu.:25.00   3rd Qu.:570.0   3rd Qu.:687.5   3rd Qu.:480.0  
##  Max.   :31.00   Max.   :740.0   Max.   :780.0   Max.   :630.0  
##  NA's   :83      NA's   :34      NA's   :34      NA's   :34     
##      VERB_3          ACT_1           ACT_3          APP_REC     
##  Min.   :440.0   Min.   :16.00   Min.   :21.00   Min.   :  787  
##  1st Qu.:520.0   1st Qu.:19.00   1st Qu.:25.00   1st Qu.: 4323  
##  Median :550.0   Median :21.00   Median :26.00   Median : 7654  
##  Mean   :561.3   Mean   :21.28   Mean   :26.55   Mean   : 8544  
##  3rd Qu.:600.0   3rd Qu.:23.00   3rd Qu.:28.00   3rd Qu.:11776  
##  Max.   :720.0   Max.   :29.00   Max.   :33.00   Max.   :48094  
##  NA's   :34      NA's   :67      NA's   :67      NA's   :1      
##     APP_ACC         NEW_STUD        NEW10           NEW25          FULLTIME    
##  Min.   :  507   Min.   : 210   Min.   : 8.00   Min.   :24.00   Min.   :  912  
##  1st Qu.: 3033   1st Qu.:1264   1st Qu.:24.00   1st Qu.:52.00   1st Qu.: 5846  
##  Median : 4761   Median :1949   Median :32.00   Median :63.00   Median :10215  
##  Mean   : 5546   Mean   :2252   Mean   :41.48   Mean   :65.55   Mean   :11296  
##  3rd Qu.: 7232   3rd Qu.:3035   3rd Qu.:57.00   3rd Qu.:82.00   3rd Qu.:15033  
##  Max.   :26330   Max.   :7425   Max.   :98.00   Max.   :99.00   Max.   :31643  
##  NA's   :1                      NA's   :16      NA's   :26                     
##     PARTTIME          IN_STATE        OUT_STAT        R_B_COST   
##  Min.   :   16.0   Min.   :  647   Min.   : 2279   Min.   :2082  
##  1st Qu.:  804.8   1st Qu.: 2100   1st Qu.: 6712   1st Qu.:3588  
##  Median : 1694.5   Median : 3030   Median : 8400   Median :4213  
##  Mean   : 2519.9   Mean   : 6641   Mean   :10108   Mean   :4538  
##  3rd Qu.: 3240.8   3rd Qu.:12348   3rd Qu.:12668   3rd Qu.:5564  
##  Max.   :21836.0   Max.   :20100   Max.   :20100   Max.   :7425  
##  NA's   :6         NA's   :7       NA's   :1                     
##       ROOM          BOARD         ADD_FEE            BOOK           PERSONAL   
##  Min.   :1033   Min.   :1000   Min.   :  20.0   Min.   : 300.0   Min.   : 300  
##  1st Qu.:1810   1st Qu.:1762   1st Qu.: 210.0   1st Qu.: 500.0   1st Qu.:1164  
##  Median :2644   Median :2125   Median : 425.5   Median : 600.0   Median :1600  
##  Mean   :2708   Mean   :2241   Mean   : 648.1   Mean   : 603.1   Mean   :1763  
##  3rd Qu.:3490   3rd Qu.:2586   3rd Qu.: 694.0   3rd Qu.: 673.8   3rd Qu.:2200  
##  Max.   :6965   Max.   :4760   Max.   :4374.0   Max.   :1230.0   Max.   :6800  
##  NA's   :42     NA's   :64     NA's   :40       NA's   :2        NA's   :11    
##       PH_D           TERM_D         SF_RATIO         DONATE     
##  Min.   :63.00   Min.   :67.00   Min.   : 2.90   Min.   : 4.00  
##  1st Qu.:80.50   1st Qu.:87.00   1st Qu.:10.88   1st Qu.:10.75  
##  Median :87.00   Median :92.00   Median :14.50   Median :17.00  
##  Mean   :85.72   Mean   :90.52   Mean   :14.23   Mean   :19.01  
##  3rd Qu.:92.00   3rd Qu.:96.00   3rd Qu.:18.02   3rd Qu.:24.00  
##  Max.   :99.00   Max.   :99.00   Max.   :24.70   Max.   :54.00  
##  NA's   :9       NA's   :16                      NA's   :12     
##     INSTRUCT        GRADUAT         SAL_FULL          SAL_AC     
##  Min.   : 3605   Min.   :10.00   Min.   : 446.0   Min.   :364.0  
##  1st Qu.: 7604   1st Qu.:47.50   1st Qu.: 597.0   1st Qu.:445.0  
##  Median : 9840   Median :62.00   Median : 661.0   Median :479.5  
##  Mean   :12832   Mean   :62.02   Mean   : 669.2   Mean   :487.1  
##  3rd Qu.:14340   3rd Qu.:74.50   3rd Qu.: 732.2   3rd Qu.:521.2  
##  Max.   :62469   Max.   :99.00   Max.   :1009.0   Max.   :733.0  
##                  NA's   :5                                       
##      SAL_AS         SAL_ALL         COMP_FUL         COMP_AC     
##  Min.   :323.0   Min.   :362.0   Min.   : 537.0   Min.   :438.0  
##  1st Qu.:381.8   1st Qu.:472.2   1st Qu.: 729.0   1st Qu.:556.5  
##  Median :407.0   Median :522.5   Median : 815.0   Median :606.0  
##  Mean   :412.8   Mean   :534.0   Mean   : 827.2   Mean   :612.1  
##  3rd Qu.:434.2   3rd Qu.:578.2   3rd Qu.: 910.0   3rd Qu.:662.2  
##  Max.   :576.0   Max.   :866.0   Max.   :1236.0   Max.   :909.0  
##                                                                  
##     COMP_AS         COMP_ALL         NUM_FULL         NUM_AC     
##  Min.   :395.0   Min.   : 436.0   Min.   : 39.0   Min.   : 32.0  
##  1st Qu.:481.0   1st Qu.: 587.0   1st Qu.:184.5   1st Qu.:138.5  
##  Median :508.5   Median : 652.0   Median :278.0   Median :208.0  
##  Mean   :519.0   Mean   : 665.8   Mean   :336.8   Mean   :230.3  
##  3rd Qu.:552.2   3rd Qu.: 730.2   3rd Qu.:457.8   3rd Qu.:299.0  
##  Max.   :717.0   Max.   :1075.0   Max.   :997.0   Max.   :721.0  
##                                                   NA's   :1      
##      NUM_AS         NUM_INS          NUM_ALL      
##  Min.   : 29.0   Min.   :  0.00   Min.   : 108.0  
##  1st Qu.:128.5   1st Qu.:  5.00   1st Qu.: 505.5  
##  Median :175.0   Median : 16.00   Median : 721.0  
##  Mean   :190.7   Mean   : 27.09   Mean   : 812.8  
##  3rd Qu.:238.2   3rd Qu.: 40.00   3rd Qu.:1035.0  
##  Max.   :510.0   Max.   :178.00   Max.   :2261.0  
##                  NA's   :1
col_I_sn$PPIND=ifelse(col_I_sn$PPIND==1, "public", "private")

NA review

col_I_sn_noNA<-col_I_sn|>na.omit()
paste("col_I_sn dim: ", paste(dim(col_I_sn),collapse=", "))
## [1] "col_I_sn dim:  176, 49"
paste("col_I_sn_noNA dim: ", paste(dim(col_I_sn_noNA), collapse=", "))
## [1] "col_I_sn_noNA dim:  22, 49"
naCol<-colSums(is.na(col_I_sn))|>sort(decreasing=TRUE)
naCol
##  AVR_ACT  AVRMATH  AVRVERB  AVRCOMB    ACT_1    ACT_3    BOARD     ROOM 
##       83       68       68       67       67       67       64       42 
##  ADD_FEE   MATH_1   MATH_3   VERB_1   VERB_3    NEW25    NEW10   TERM_D 
##       40       34       34       34       34       26       16       16 
##   DONATE PERSONAL     PH_D IN_STATE PARTTIME  GRADUAT     BOOK  APP_REC 
##       12       11        9        7        6        5        2        1 
##  APP_ACC OUT_STAT   NUM_AC  NUM_INS     ...1    PPIND     FICE    STATE 
##        1        1        1        1        0        0        0        0 
##     TYPE NEW_STUD FULLTIME R_B_COST SF_RATIO INSTRUCT SAL_FULL   SAL_AC 
##        0        0        0        0        0        0        0        0 
##   SAL_AS  SAL_ALL COMP_FUL  COMP_AC  COMP_AS COMP_ALL NUM_FULL   NUM_AS 
##        0        0        0        0        0        0        0        0 
##  NUM_ALL 
##        0
table(naCol)
## naCol
##  0  1  2  5  6  7  9 11 12 16 26 34 40 42 64 67 68 83 
## 21  5  1  1  1  1  1  1  1  2  1  4  1  1  1  3  2  1

log dataset

q_cols<-sapply(col_I_sn, is.numeric)
q_cols[c("PPIND", "FICE")]<-FALSE
q_cols<-names(q_cols)[q_cols==TRUE]

col_I_sn_log<-col_I_sn|>
  mutate(across(q_cols, ~log(.x)))|>
  rename_with(~paste0("log-", .x), q_cols)

head(col_I_sn_log)|>rmarkdown::paged_table()

1.1 Разобраться в том, что означают признаки.

Выделим основные признаки, по которым будет проходить дальнейший анализ данных.

1.1.1 Список признаков.

1.1.1.1 Количественные признаки:

  • AVRMATH Average Math SAT score
  • AVRVERB Average Verbal SAT score
  • AVRCOMB Average Combined SAT score
  • AVR_ACT Average ACT score
  • MATH_1 First quartile - Math SAT
  • MATH_3 Third quartile - Math SAT
  • VERB_1 First quartile - Verbal SAT
  • VERB_3 Third quartile - Verbal SAT
  • ACT_1 First quartile - ACT
  • ACT_3 Third quartile - ACT
  • APP_REC Number of applications received
  • APP_ACC Number of applicants accepted
  • NEW_STUD Number of new students enrolled
  • FULLTIME Number of fulltime undergraduates
  • PARTTIME Number of parttime undergraduates
  • IN_STATE In-state tuition
  • OUT_STAT Out-of-state tuition
  • R_B_COST Room and board costs
  • ROOM Room costs
  • BOARD Board costs
  • ADD_FEE Additional fees
  • BOOK Estimated book costs
  • PERSONAL Estimated personal spending
  • PH_D Pct. of faculty with Ph.D.’s
  • TERM_D Pct. of faculty with terminal degree
  • SAL_FULL Average salary - full professor
  • SAL_AC Average salary - associate professor
  • SAL_AS Average salary - assistant professor
  • SAL_ALL Average salary - all ranks
  • COMP_FUL Average compensation - full professor
  • COMP_AC Average compensation - associate professor
  • COMP_AS Average compensation - assistant professor
  • COMP_ALL Average compensation - all ranks
  • NUM_FULL Number of full professor
  • NUM_AC Number of associate professor
  • NUM_AS Number of assistant professor
  • NUM_INS Number of instructors
  • NUM_ALL Number of faculty - all ranks
  • INSTRUCT Instructional expenditure per student
  • GRADUAT Graduation rate
  • SF_RATIO Student/faculty ratio
  • DONATE Pct.alumni who donate
  • NEW10 Pct. new students from top 10% of H.S. class - % студентов из топ 10% своей старшей школы
  • NEW25 Pct. new students from top 25% of H.S. class - % студентов из топ 25% своей старшей школы

1.1.1.2 Качественные признаки:

  • FICE - Federal ID Number
  • …1 - Название университета
  • PPIND Public/private indicator (public=1, private=2)
  • STATE State (postal code)
  • TYPE - I (можно удалить)

Всего 49 признаков.

1.2 Если признаков очень много, то отобрать признаки (примерно 7-10)…

Отберём из условий задачи и по смыслу

  1. PPIND (для группировки)

Из условий задачи: 2. ADD_FEE - дополнительные выплаты 3. BOOK - плата за бронирование 4. NEW10 (зависимая переменная, также NEW25)

NEW10 зависит от AVRMATH, MATH_1, MATH_3, AVRVERB, VERB_1, VERB_3, AVR_ACT, ACT_1, ACT_3 и AVRCOMB. Признаков много, так что выберем самый обобщающий - AVRCOMB 5. AVRCOMB

  1. SF_RATIO

  2. PH_D

  3. GRADUAT 6-8 из предположения, что они взаимосвязаны

  4. R_B_COST

  5. IN_STATE

  6. OUT_STAT 9-11 - зависимость стоимости проживания (как окажется дальше из matrix plot, IN_STATE - фиктивный признак, то есть не информативный)

1.3 Определить вид признаков (колич., порядковые, качеств.)…

Проделаем это для рассматриваемых 11 признаков.

Setup

colNames<-c("ADD_FEE", "BOOK", "R_B_COST","IN_STATE","OUT_STAT","NEW10","AVRCOMB","SF_RATIO","PH_D","GRADUAT","PPIND")
colFin<-c("ADD_FEE", "BOOK", "PPIND")
colCost<-c("R_B_COST", "IN_STATE", "OUT_STAT", "PPIND")
colNew<-c("NEW10", "AVRCOMB", "PPIND")
colStud<-c("SF_RATIO", "PH_D", "GRADUAT", "PPIND")
col_I_comp<-col_I_sn[colNames]
col_I_comp_log<-col_I_sn_log[c(paste0("log-", colNames)[-length(colNames)], "PPIND")]
head(col_I_comp_log)|>rmarkdown::paged_table()
unique_info <- function(column) {
  column<-na.omit(column)
  list(ratio=length(unique(column)) / length(column), unique=length(unique(column)), total=length(column), moda = max(table(column)))
}

ratios <- sapply(col_I_comp, unique_info)
ratios
##        ADD_FEE BOOK      R_B_COST  IN_STATE  OUT_STAT  NEW10 AVRCOMB  
## ratio  0.875   0.3793103 0.9431818 0.9053254 0.9314286 0.425 0.7889908
## unique 119     66        166       153       163       68    86       
## total  136     174       176       169       175       160   109      
## moda   4       29        3         4         4         10    3        
##        SF_RATIO  PH_D      GRADUAT   PPIND     
## ratio  0.6363636 0.1976048 0.3918129 0.01136364
## unique 112       33        67        2         
## total  176       167       171       176       
## moda   5         10        7         123

1.3.1 Количественные признаки

1.3.1.1 Непрерывные признаки

  1. ADD_FEE
  2. R_B_COST
  3. IN_STATE
  4. OUT_STATE
  5. AVRCOMB
  6. SF_RATIO
  7. GRADUAT

1.3.1.2 Дискретные признаки

  1. BOOK (непрерывный, но дискретный из-за моды. Может быть округление)
  2. PH_D (считаем из-за плохой точности и округления данных)
  3. NEW10 (по ratio и moda, низкая точность)

1.3.2 Качественные признаки

  • PPIND - private/public.

1.4 Построить matrix plot (pairs plot), его долго разглядывать с точки зрения outliers, неоднородностей, вида распределений, вида зависимостей (линейные/нелинейные) и пр.

1.4.1 Первое построение pairs plot. Общие наблюдения.

ggpairs(col_I_comp, aes(alpha = 0.5, color=PPIND),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))

1.4.1.1 Наблюдения:

  1. IN_STATE и OUT_STAT для private выстраивает линейную зависимость, для public - близкая к линейной зависимости. Отличительная особенность признака IN_STATE от OUT_STAT - наличие других факторов на образование цены обучения. OUT_STATE более приложим для анализа данных, поэтому IN_STATE можно считать фиктивным признаком.
  2. ADD_FEE похоже на log-normal. Заметны Outliers.
1.4.1.1.1 Outliers, ADD_FEE:
col_I_sn|>filter(ADD_FEE>3000)

University of California, Massachusetts - особые индивиды, в связи с высоким ВВП штатов.

  1. По histogram plot неоднородности (по PPIND) в R_B_COST, IN_STATE, OUT_STAT, SF_RATIO, GRADUAT.

Видим, что NEW10 зависит от PH_D, GRADUAT, AVRCOMB, отдельно для private зависит от OUT_STAT, IN_STATE, для public - от ADD_FEE, BOOK, R_B_COST

1.5 Если есть сильно несимметричные (с хвостом вправо) распределения на положительной полуоси, то…

1.5.1 Второе построение pairs plot.

Прологарифмируем ADD_FEE

col_I_mixed<-col_I_comp|>mutate(ADD_FEE=col_I_comp_log$`log-ADD_FEE`)|>rename(`log-ADD_FEE`=ADD_FEE)|>dplyr::select(-IN_STATE)
head(col_I_mixed)
ggpairs(col_I_mixed, aes(alpha = 0.5, color=PPIND),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

1.5.1.1 Наблюдения:

  • Неоднородности (по PPIND) R_B_COST, OUT_STAT, NEW10.
  • Близкие к линейномй зависимости NEW10/ AVRCOMB.
  • Outliers, BOOK
1.5.1.1.1 Outliers, BOOK

Особые индивиды

col_I_sn|>filter(BOOK>=900 | BOOK <= 350)

Howard, Tulane, Hofstra University, Rensselaer Polytechn - особенности расположения, штата.

1.6 Если есть outliers, то попробовать объяснить причину (ошибка в данных, особые индивиды) и удалить их.

1.6.1 BOOK

ggpairs(col_I_mixed[,c("log-ADD_FEE", "BOOK", "PPIND")], aes(alpha = 0.5, color=PPIND),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

col_I_sn|>filter(BOOK>=900 | BOOK <= 350)

Может быть как низкой, так и высокой в зависимости от штата, политики университета и скрытых факторов.

1.6.2 ADD_FEE

col_I_sn|>filter(ADD_FEE>3000)

University of California, Massachusetts - особые индивиды (из-за ВВП).

1.6.3 AVRCOMB/NEW10

ggpairs(col_I_mixed[colNew], aes(alpha = 0.5, color=PPIND),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

col_I_sn|>filter(AVRCOMB<=1200 & NEW10>=75)

University of California, North Carolina. Причиной может быть низкий порог для поступления.

1.6.4 PH_D

ggpairs(col_I_mixed[colStud], aes(alpha = 0.5, color=PPIND),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

col_I_sn|>filter(PH_D<=70)

Незначительные выбросы - особые индивиды. Могут быть связаны с низким GRADUAT и высоким SF_RATIO.

1.6.5 OUT_STAT

ggpairs(col_I_mixed[,c("R_B_COST", "OUT_STAT", "PPIND")], aes(alpha = 0.5, color=PPIND),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

col_I_sn|>filter(OUT_STAT>=14000 & PPIND=="public")
col_I_sn|>filter(OUT_STAT<=10500 & PPIND=="private")

Цены на обучение варьируются по штатам и самим университетам.


1.7 Если есть неоднородности (например, видны два облака точек), то объяснить причину (найти категоризующую переменную, объясняющую эту неоднородность).

ggpairs(col_I_mixed, aes(alpha = 0.5, color=PPIND),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Выделим неоднородности R_B_COST, OUT_STAT. Явное разделение на два облака точек для этих признаков может быть связано с разницей финансирования. Для public - государство, штат и граждане, для private - только граждане.


1.8 В дальнейшем вид pairs plots, распределения признаков и корреляции анализировать отдельно для неоднородных групп.


1.9 Посмотрите также на descriptive statistics с точки зрения минимумов-максимумов, асимметрии, эксцесса и пр.

results <- describeBy(col_I_mixed, group=col_I_comp$PPIND)
df_results <- map_dfr(results, ~as.data.frame(.x), .id = "group")
kable(df_results, "html") %>%
  kable_styling(full_width = FALSE) %>%
  scroll_box(width = "100%", height = "400px")
group vars n mean sd median trimmed mad min max range skew kurtosis se
log-ADD_FEE…1 private 1 46 5.838785 0.7936613 6.003811 5.872662 0.7083206 3.688880 7.515345 3.826465 -0.4493357 0.0057261 0.1170190
BOOK…2 private 2 53 610.075472 158.1109456 600.000000 597.534884 148.2600000 300.000000 1230.000000 930.000000 1.2325851 3.0965764 21.7182087
R_B_COST…3 private 3 53 5839.811321 997.0704793 5975.000000 5922.255814 816.9126000 3320.000000 7425.000000 4105.000000 -0.7233444 -0.2295138 136.9581633
OUT_STAT…4 private 4 53 15747.358491 3924.9053349 17020.000000 16306.186047 2772.4620000 2340.000000 20100.000000 17760.000000 -1.3082365 1.3192370 539.1272103
NEW10…5 private 5 52 57.788461 24.3796255 58.000000 57.761905 33.3585000 16.000000 98.000000 82.000000 0.0381849 -1.4295240 3.3808458
AVRCOMB…6 private 6 29 1177.034483 144.2860459 1199.000000 1177.680000 161.6034000 883.000000 1410.000000 527.000000 -0.0893581 -1.2296542 26.7932461
SF_RATIO…7 private 7 53 9.813208 4.3268854 9.200000 9.555814 4.8925800 2.900000 20.500000 17.600000 0.4758614 -0.7282486 0.5943434
PH_D…8 private 8 49 89.591837 7.8658502 91.000000 90.317073 7.4130000 71.000000 99.000000 28.000000 -0.8078233 -0.4898133 1.1236929
GRADUAT…9 private 9 52 78.711539 15.7534142 78.500000 80.023809 17.7912000 33.000000 99.000000 66.000000 -0.6127581 -0.1790876 2.1846055
PPIND…10 private 10 53 1.000000 0.0000000 1.000000 1.000000 0.0000000 1.000000 1.000000 0.000000 NaN NaN 0.0000000
log-ADD_FEE…11 public 1 90 6.014273 1.1296149 6.093567 6.024473 0.8723448 2.995732 8.383433 5.387701 -0.1980167 0.2896088 0.1190719
BOOK…12 public 2 121 600.008264 101.1621220 600.000000 591.103093 111.1950000 400.000000 858.000000 458.000000 0.6447669 -0.2889320 9.1965565
R_B_COST…13 public 3 123 3977.073171 842.4904744 3811.000000 3911.979798 722.0262000 2082.000000 6607.000000 4525.000000 0.7886999 0.7902446 75.9648078
OUT_STAT…14 public 4 122 7657.549180 2252.3124723 7446.500000 7517.163265 1939.2408000 2279.000000 15732.000000 13453.000000 0.7771996 1.2996173 203.9147900
NEW10…15 public 5 108 33.620370 21.8669373 26.000000 29.920455 11.8608000 8.000000 95.000000 87.000000 1.5420644 1.5799521 2.1041470
AVRCOMB…16 public 6 80 1014.662500 84.6586510 997.000000 1010.687500 75.6126000 810.000000 1240.000000 430.000000 0.4180047 0.1997793 9.4651249
SF_RATIO…17 public 7 123 16.129268 3.8302021 16.500000 16.178788 4.1512800 6.700000 24.700000 18.000000 -0.0994206 -0.5699091 0.3453577
PH_D…18 public 8 118 84.110169 7.7985139 85.000000 84.406250 7.4130000 63.000000 99.000000 36.000000 -0.3698537 -0.3490528 0.7179114
GRADUAT…19 public 9 119 54.731092 15.5730632 54.000000 54.412371 16.3086000 10.000000 95.000000 85.000000 0.1127780 -0.1892703 1.4275804
PPIND…20 public 10 123 2.000000 0.0000000 2.000000 2.000000 0.0000000 2.000000 2.000000 0.000000 NaN NaN 0.0000000

1.9.1 Наблюдения:

1.9.1.1 Mean, median, skew:

NEW10 и AVRCOMB имеют skew, близкий к 0 (private), SF_RATIO, GRADUAT, log-ADD_FEE околонулевой skew (public).

1.9.1.2 Kurtosis:

Близкий к 0 (private): log-ADD_FEE (но skew ~ -0.44).

2 2. О виде распределений и о сравнении распределений

2.1 Первые два пункта индивидуального задания нужно делать не по указанному порядку, а как того требует логика статистики. Чтобы сравнивать выборки по t-критерию, нужно знать о том, близки ли распределения в сравниваемых группах к нормальным или хотя бы к унимодальным и симметричным. Чтобы проверять распределения признаков на нормальность, нужно знать, что рассматривается однородная выборка.

col_split_mixed<-split(col_I_mixed, col_I_mixed$PPIND)
names(col_split_mixed)<-c("private", "public")
colMixed<-names(col_I_mixed)
colMixed<-colMixed[-length(colMixed)]
colMixed
## [1] "log-ADD_FEE" "BOOK"        "R_B_COST"    "OUT_STAT"    "NEW10"      
## [6] "AVRCOMB"     "SF_RATIO"    "PH_D"        "GRADUAT"

2.1.1 Проверка на нормальность по хи-квадрат критерию (Авторство: Яковлев Д.М.)

chi_squared_normality_test <- function(data) {
  data <- na.omit(data)
  n <- length(data)
  if (n < 10) stop("Sample size too small for this test. Need at least 10 observations.")

  mu <- mean(data)
  sigma <- sd(data)

   # Calculate quantiles for breaks
   probs <- seq(5/n, 1 - 5/n, by = 5/n)
   breaks <- qnorm(probs, mean = mu, sd = sigma)
  
  breaks <- c(qnorm(0), breaks, qnorm(1))
  #Observed Frequencies
  observed <- hist(data, breaks = breaks, plot = FALSE)$counts

  #Expected Frequencies (using pnorm)
  expected <- diff(pnorm(breaks, mean = mu, sd = sigma)) * n
    
  chi_squared <- sum((observed - expected)^2 / expected)
  df <- length(expected) - 3  #Corrected degrees of freedom

  p_value <- pchisq(chi_squared, df = df, lower.tail = FALSE)

  return(list(p_value = p_value, statistic = chi_squared, df = df))
}

2.1.1.1 Проверка на достоверность (p-value должно иметь равномерное распределение на (0, 1))

n = 500
test<-replicate(5000, chi_squared_normality_test(rnorm(n))$p_value)
quantile(test)
##          0%         25%         50%         75%        100% 
## 0.000116153 0.249925966 0.509619748 0.761571168 0.999679814

2.1.2 Тесты на нормальность

perform_normality_tests <- function(df, cols) {
  # Input validation (remains the same)
    # Function to perform a single normality test (simplified)
  run_test <- function(data, test_name) {
    if (test_name == "Chi-squared") {
      chisq_test<-chi_squared_normality_test(data)
      return(chisq_test[c("p_value", "statistic", "df")])
    }
    test_result <- switch(test_name,
                          "Lilliefors" = lillie.test(data),
                          "Anderson-Darling" = ad.test(data),
                          "Shapiro-Wilk" = shapiro.test(data)
    )
    return(list(p_value = test_result$p.value, statistic = test_result$statistic, df=length(na.omit(data))))  # Removed test_name
  }
     results <- lapply(cols, function(col) {
    data <- df[[col]]
    n_na <- sum(is.na(data))
    if (n_na > 0) {
      warning(paste0("Removed ", n_na, " NA values from column '", col, "' before testing."))
      data <- na.omit(data)
    }
    test_results <- lapply(c("Lilliefors", "Anderson-Darling", "Shapiro-Wilk", "Chi-squared"), function(test) run_test(data, test))
    list(column = col, tests = test_results)
  })

  results_df <- do.call(rbind, lapply(results, function(x) {
    test_names <- c("Lilliefors", "Anderson-Darling", "Shapiro-Wilk", "Chi-squared")
    p_values <- sapply(x$tests, "[[", "p_value")
    data.frame(
      Column = x$column,
      Test = test_names, #Directly use test_names vector
      Statistic = sapply(x$tests, "[[", "statistic"),
      P_value = p_values,
      Size = sapply(x$tests, "[[", "df"),
      Significance = sapply(p_values, function(p) ifelse(p < 0.05, "\u2714", "\u2716")),
      stringsAsFactors = FALSE
    )
  }))

  plots <- lapply(cols, function(col) {
    data <- df[[col]]
    data <- na.omit(data)
    ggqqplot(data, title = paste("Normal Probability Plot for", col))

  })
  names(plots)<-cols
  list(results = results_df, plots = plots)
}

2.2 Так как визуально однородность при предварительном анализе была уже исследована, то можно начинать с анализа вида распределения признаков, возможно, по группам. Сюда входит: normal probability plot (что это такое?), проверка по критериям Лиллиефорса, AD, хи-квадрат, Шапиро-Уилка. По критерию хи-квадрат, а также визуально по PP-plot можно проверить и гипотезы о согласии с другими распределениями, например, логнормальным. Задаваемые вопросы: чем отличается критерий Лиллиефорса от критерия Колмогорова, в чем отличие AD, как выглядят статистики критериев, что такое PP-plot и normal probability plot, почему естественно при рисовании normal probability plot одновременно смотреть на результаты критерия Шапиро-Уилка.

2.2.1 Public: проверка на нормальность, тесты

# colSums(is.na(col_split_mixed$public))
col_norm_public<-col_split_mixed$public|>perform_normality_tests(colMixed)
col_norm_public$plots
## $`log-ADD_FEE`

## 
## $BOOK

## 
## $R_B_COST

## 
## $OUT_STAT

## 
## $NEW10

## 
## $AVRCOMB

## 
## $SF_RATIO

## 
## $PH_D

## 
## $GRADUAT

col_norm_public$results|>dplyr::select(-Statistic)|>rmarkdown::paged_table()

2.2.1.1 Наблюдения:

Хи-квадрат: не отвергаем AVRCOMB, SF_RATIO, PH_D, GRADUAT, R_B_COST, OUT_STAT,

Можно заметить “хвосты” на некоторых графиках.

2.2.2 Private: проверка на нормальность, тесты

col_norm_private<-col_split_mixed$private|>perform_normality_tests(colMixed)
col_norm_private$plots
## $`log-ADD_FEE`

## 
## $BOOK

## 
## $R_B_COST

## 
## $OUT_STAT

## 
## $NEW10

## 
## $AVRCOMB

## 
## $SF_RATIO

## 
## $PH_D

## 
## $GRADUAT

col_norm_private$results|>dplyr::select(-Statistic)|>rmarkdown::paged_table()

2.2.2.1 Наблюдения:

Хи-квадрат: не отвергаем log-ADD_FEE, R_B_COST, NEW10, AVRCOMB, SF_RATIO, GRADUAT

2.3 Сначала имеет смысл посмотреть на сравнение распределений в группах с помощью ящиков с усами. С помощью ящиков с усами там, где групп больше двух, можно выбрать две из них, которые интересно сравнить с помощью критериев.

Совместно со следующим заданием.


2.4 Если в задании есть сравнение независимых выборок (точнее, распределений независимых случайных величин), то начинать нужно с t-критерия, который мощный против альтернатив, заключающихся в наиболее легко интерпретируемом сдвиге (т.е. разнице средних). Нужно не забыть, что у критерия есть варианты для модели с одинаковыми дисперсиями (получается точное p-value, которое может быть неправильным, если на самом деле дисперсии одинаковые) и с произвольными дисперсиями. В результате, получатся два критерия для гипотезы о равенстве средних и два критерия о равенстве разбросов. Нужно уметь объяснять, что это за критерии и при каких условиях их можно применять. Не забудьте, что при использовании асимптотических критериев нужно обращать внимание на объемы выборок. Сделайте выводы о том, для каких признаков есть разница в сдвиге.

2.4.1 Тесты на равенство дисперсий/матожиданий, тест Wilcoxon.

perform_tests <- function(data,
                          quantitative_columns,
                          categorical_column,
                          alpha = 0.05) {
  results <- list()
  
  for (col in quantitative_columns) {
    levene_test<-leveneTest(data[[col]] ~ data[[categorical_column]])
    var_test <- var.test(data[[col]] ~ data[[categorical_column]])
    wilcox_test <- wilcox.test(data[[col]]~data[[categorical_column]])
    t_test_equal <- t.test(data[[col]] ~ data[[categorical_column]], var.equal = TRUE)
    t_test_unequal <- t.test(data[[col]] ~ data[[categorical_column]], var.equal = FALSE)
    ks_test<-ks.test(data[[col]]~data[[categorical_column]])
  
    results[[col]] <- list(
      t_test_unequal_p_value = t_test_unequal$p.value,
      t_test_equal_p_value = t_test_equal$p.value,
      wilcox_test_p_value = wilcox_test$p.value,
      levene_test_p_value = levene_test$`Pr(>F)`[1],
      var_test_p_value = var_test$p.value,
      ks_test_p_value = ks_test$p.value
    )
  }
  
  return(results |> enframe(name = "attribute", value = "p_values") |>
  unnest_wider(p_values))
}
df_tests = perform_tests(col_I_mixed, colMixed, "PPIND")
df_long <- col_I_mixed|>
  pivot_longer(cols = -PPIND, names_to = "attribute", values_to = "value")
# Объединяем результаты теста с основными данными
df_with_tests <- df_long %>%
  left_join(df_tests, by = "attribute")
df_with_tests_cost <- df_with_tests[df_with_tests$attribute%in%colCost[-length(colCost)],]
df_with_tests_new <- df_with_tests[df_with_tests$attribute%in%colNew[-length(colNew)],]
df_with_tests_stud <- df_with_tests[df_with_tests$attribute%in%colStud[-length(colStud)],]
df_with_tests_fin <- df_with_tests[df_with_tests$attribute%in%c("log-ADD_FEE", "BOOK"),]

p<-ggplot(df_with_tests_cost,
       aes(x = attribute,
           y = value,
           fill = PPIND)) +
  geom_boxplot() +
  labs(title = "Boxplot",
       x = "Attribute",
       y = "Value")+
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), plot.caption = element_text(size = 8, hjust = 0)) +
  scale_fill_manual(values = c("private" = "red", "public" = "lightblue"))

2.4.2 AVRCOMB, NEW10.

p$data<-df_with_tests_new|>filter(attribute%in%"AVRCOMB")
p

p$data<-df_with_tests_new|>filter(attribute%in%"NEW10")
p

df_with_tests|>dplyr::filter(attribute%in%colNew[-length(colNew)])|>dplyr::select(-PPIND, -value)|>slice_head(n=2)|>rmarkdown::paged_table()

2.4.2.1 Наблюдения:

Можем сказать, что public и private для AVRCOMB и NEW10 имеют разные матожидания.

В случае гипотезы о равенстве дисперсий, критерий Левена (Levene) и критерий Фишера (var) имеют разные мощности. Критерий Фишера используется, если сравниваемые выборки имеют нормальное распределение. Поскольку группа public признака NEW10 отвергла гипотезу о нормальном распределении, применять критерий Фишера нельзя.

Отвергаем равенство дисперсий для AVRCOMB, NEW10.

2.4.3 BOOK, log-ADD_FEE.

p$data<-df_with_tests_fin|>filter(attribute%in%"BOOK")
p

p$data<-df_with_tests_fin|>filter(attribute%in%"log-ADD_FEE")
p

df_with_tests|>dplyr::filter(attribute%in%c("log-ADD_FEE", "BOOK"))|>dplyr::select(-PPIND, -value)|>head(length(colFin)-1)|>rmarkdown::paged_table()

2.4.3.1 Наблюдения:

Исходя из p_value не отвергаем, что public и private имеют равные матожидания. Тест Вилкоксона также характеризует их однородность.

2.4.3.2 Задача 1.1 из .tsk

Проверим log-ADD_FEE и BOOK друг с другом

compare_column_tests <- function(column_1,
                          column_2,
                          alpha = 0.05) {
  levene_test<-leveneTest(column_2, column_1)
  var_test <- var.test(column_1, column_2)
  wilcox_test <- wilcox.test(column_1, column_2)
  t_test_equal <- t.test(column_1, column_2, var.equal = TRUE)
  t_test_unequal <- t.test(column_1, column_2, var.equal = FALSE)
  
  results <- data.frame(
    t_test_unequal_p_value = t_test_unequal$p.value,
    t_test_equal_p_value = t_test_equal$p.value,
    levene_test_p_value = levene_test$`Pr(>F)`[1],
    var_test_p_value = var_test$p.value,
    wilcox_test_p_value = wilcox_test$p.value
  )
}
af_pr<-col_I_sn$ADD_FEE[col_I_sn$PPIND=="private"]
af_pu<-col_I_sn$ADD_FEE[col_I_sn$PPIND=="public"]
b_pr<-col_split_mixed$private$BOOK
b_pu<-col_split_mixed$public$BOOK

print(compare_column_tests(af_pr, b_pr))
##   t_test_unequal_p_value t_test_equal_p_value levene_test_p_value
## 1            0.005903872          0.003538701        2.925758e-75
##   var_test_p_value wilcox_test_p_value
## 1      1.42501e-07        1.678377e-05
print(compare_column_tests(af_pu, b_pu))
##   t_test_unequal_p_value t_test_equal_p_value levene_test_p_value
## 1              0.1706573            0.1115642           0.9495775
##   var_test_p_value wilcox_test_p_value
## 1                0        0.0006091651
2.4.3.2.1 Наблюдения:

В общем случае тесты показывают, что разница есть в случае общественных вузов, но для частных вузов не отвергаем равенство дисперсий (из levene’s test) и матожиданий.

2.4.4 GRADUAT, PH_D, SF_RATIO

p$data<-df_with_tests_stud
p

df_with_tests|>dplyr::filter(attribute%in%colStud[-length(colStud)])|>dplyr::select(-PPIND, -value)|>head(length(colStud)-1)|>rmarkdown::paged_table()

2.4.4.1 Наблюдения:

Отвергаем равенство матожиданий для GRADUAT, SF_RATIO, PH_D. Не отвергаем равенство дисперсий GRADUAT, SF_RATIO, PH_D. — ### OUT_STAT, R_B_COST

p$data<-df_with_tests_cost
p

df_with_tests|>dplyr::filter(attribute%in%colCost[-length(colCost)])|>dplyr::select(-PPIND, -value)|>head(length(colCost)-2)|>rmarkdown::paged_table()

2.4.4.2 Наблюдения:

Отвергаем равенство матожиданий для OUT_STAT, R_B_COST. Не отвергаем равенство дисперсий для R_B_COST.

2.5 Далее, объясняете, в каких случаях (распределение далеко от нормального, могут быть выделяющиеся наблюдения, нас интересует другая характеристика положения) t-критерий не удовлетворителен и нужно переходить к непараметрическим критериям. Рассказываете, какой из непараметрических критериев является аналогом t-критерия, как он строится и против какой альтернативы мощный. Вы уже догадались, что это критерий Манна-Уитни, он же критерий Вилкоксона.

t-критерий неудовлетворителен в нескольких ситуациях:

  1. Малый объём выборки. Если форма распределения признака неизвестна или не нормальна, то t-критерий асимптотически точный.
  2. Наличие выбросов. t-критерий неусточив по отношению к выбросам.
  3. Другие цели эксперимента. t-критерий мощный против альтернатив о неравенстве матожиданий, его нельзя применить для построения других гипотез, например, гипотезы о равенстве дисперсий.

В качестве замены t-критерия можно воспользоваться непараметрическим критерием Манна-Уитни (Вилкоксона).

2.6 Смотрите на результаты применения критерия Манна-Уитни, сравниваете с результатами применения t-критерия. Проводите сравнительный анализ критериев с теоретической точки зрения (чем один лучше другого и чем хуже).

df_tests|>dplyr::select(attribute,t_test_unequal_p_value, t_test_equal_p_value, wilcox_test_p_value)|>rmarkdown::paged_table()

2.6.1 Наблюдения:

В случаях, где Манн-Уитни мощнее t-test: BOOK, R_B_COST, OUT_STAT, SF_RATIO, GRADUAT; причиной этого могут быть выбросы. В случаях, где Манн-Уитни не мощнее t-test: log-ADD_FEE, NEW10, AVRCOMB, PH_D; могут отсутствовать выбросы, из-за чего t-test оказывается мощнее.

2.6.2 Теоретическое сравнение:

Как непараметрический критерий, тест Вилкоксона уступает t-тесту в мощности против альтернативы о неравенстве матожиданий. С другой стороны, тест Вилкоксона устойчив к выбросам, мощный против альтернативы о неравенстве распределения двух случайных величин. Если известно, что две случайные величины имеют одинаковые формы распределения, то тест Вилкоксона автоматически ставит гипотезу о равенстве матожиданий.

2.7 Далее, переходите к критериям, которые умеют сравнивать не только характеристики положения, но и формы распределений (это критерий Колмогорова-Смирнова, например). Для каждого критерия (включая критерий Манна-Уитни), нужно уметь объяснять, что означают столбцы в таблицах результатов критериев. Также, при разных результатах проверки гипотезы о равенстве распределений нужно объяснять, почему один критерий, например, не отверг гипотезу, а другой – отверг.

2.7.1 Сравнение форм распределения public и private

df_tests|>dplyr::select(attribute, wilcox_test_p_value, ks_test_p_value)|>rmarkdown::paged_table()

2.7.1.1 Наблюдения:

Столбец wilcox_test_p_value - p-value из критерия Вилкоксона. Столбец ks_test_p_value - p-value из критерия Колмогорова-Смирнова. Критерии Вилкоксона и Колмогорова-Смирнова мощны против альтернативы о неравенстве распределения. Критерии не отвергают равенство форм распределений для log-ADD_FEE и BOOK, но отвергают для остальных. Среди рассматриваемых признаков нет случаев, когда один критерий отвергает гипотезу, а другой - не отвергает.


3 Об анализе зависимостей

3.1 Вспомните, какие бывают виды зависимостей и чем они измеряются, по каким формулам. Посмотрите на основе pairs plot, какие зависимости у вас в данных. Не забудьте, что при неоднородных данных изучать зависимости имеет смысл только внутри групп по-отдельности.

col_split_mixed<-split(col_I_mixed, col_I_mixed$PPIND)
names(col_split_mixed)<-c("private", "public")
col_split_mixed<-lapply(col_split_mixed, function(x){
  dplyr::select(x, -PPIND)
})
colMixed<-names(col_I_mixed)
colMixed<-colMixed[-length(colMixed)]
col_split_mixed$private|>rmarkdown::paged_table()
col_split_mixed$public|>rmarkdown::paged_table()

3.1.1 Private:

ggpairs(col_split_mixed$private, aes(alpha = 0.5),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))

3.1.1.1 Наблюдения:

Можно заметить, что NEW10/AVRCOMB разбивается на две части (поскольку AVRCOMB имеет много пропусков и сам AVRCOMB подразумевает разделение на подгруппы).

col_I_sn[col_I_sn$PPIND=="private",]|>filter(NEW10<=50 & !is.na(AVRCOMB))|>rmarkdown::paged_table()
col_I_sn[col_I_sn$PPIND=="private",]|>filter(NEW10>50 & !is.na(AVRCOMB))|>rmarkdown::paged_table()

Среди частных университетов только один “Delaware University” имеет расхождения между IN_STATE и OUT_STAT. Сейчас идёт такое ценообразование.

3.1.1.1.1 University of Delaware.

University of Delaware’s tuition is $15,410 for in-state and $37,930 for out-of-state.

col_I_sn[col_I_sn$IN_STATE!=col_I_sn$OUT_STAT & col_I_sn$PPIND=="private",]

3.1.2 Public:

ggpairs(col_split_mixed$public, aes(alpha = 0.5),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))

3.1.2.1 Наблюдения

На pairs plot можно заметить выбросы в признаках log-ADD_FEE, OUT_STAT, NEW10, которые влияют на коэффициент корреляции Пирсона.

3.2 Начинать нужно с анализа линейных зависимостей. На основе коэффициента корреляции Пирсона нужно проинтерпретировать значимые зависимости. При наличие в данных пропусков обратите внимание на выбор между casewise and pairwise MD deletion (в чем разница, какие недостатки и достоинства у этих вариантов?).

Учитывая ограниченный объём данных, тем более после разделения по группам (например, наименьшая в private размера 29 без NA), следует воспользоваться pairwise. Более того, из-за признака AVRCOMB использование casewise неявно рассматривает корреляцию в другой подгруппе. Минусом pairwise является то, что из-за попарной проверки могут учитываться данные, пропущенные в других парах.


3.3 Затем можно переходить к ранговым коэффициентам корреляции. Расскажите, при каких условиях коэффициенты корреляции Пирсона и Спирмена примерно равны. Приведите примеры, когда один из них больше другого и наоборот. Сравните результаты на ваших данных. Если при сравнении будут найдены заметные различия в результатах, то попробуйте объяснить причину.

compute_corr <- function(group_data) {
  corr_pairwise_pearson <- cor(group_data, use = "pairwise.complete.obs", method = "pearson")
  corr_pairwise_spearman <- cor(group_data, use = "pairwise.complete.obs", method = "spearman")
  return(list(pairwise_pearson = corr_pairwise_pearson, pairwise_spearman = corr_pairwise_spearman))
}

make_corr_plots <- function(corr_list, group) {
  p1 <- ggcorrplot(corr_list$pairwise_pearson,
                   lab = TRUE,
                   lab_size = 1.9,
                   colors = c("blue", "white", "red"),
                   title = paste("Pairwise Cor Pearson Matrix for Group", group),
                   ggtheme = theme_minimal() +
                     theme(plot.title = element_text(size = 6)) # Reduced title size
  )
  p2 <- ggcorrplot(corr_list$pairwise_spearman,
                   lab = TRUE,
                   lab_size = 1.9,
                   colors = c("blue", "white", "red"),
                   title = paste("Pairwise Cor Spearman Matrix for Group", group),
                   ggtheme = theme_minimal() +
                     theme(plot.title = element_text(size = 6)) # Reduced title size
  )
  grid.arrange(p1, p2, ncol = 2, widths = c(4, 4),  # Adjust widths here
               padding = unit(1, "cm"))  # Add padding for extra space
}
library(ggcorrplot)
for (group in names(col_split_mixed)) {
    group_data <- col_split_mixed[[group]]
    cor_matrix <- compute_corr(group_data)
    plot_pearson_spearman <- make_corr_plots(cor_matrix, group)
    print(plot_pearson_spearman)
}

## TableGrob (1 x 2) "arrange": 2 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]

## TableGrob (1 x 2) "arrange": 2 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]

3.3.1 Наблюдения, различия между коэффициентом корреляции Пирсона и Спирмена:

3.3.1.1 Private:

Пирсон мощнее Спирмена:

  • OUT_STAT/R_B_COST - почти линейный график и выбросы,
  • SF_RATIO/BOOK - выбросы при данных, похожих на облако.
  • R_B_COST/PH_D - выбросы усиливают линейную зависимость,
  • R_B_COST/GRADUAT - выбросы усиливают линейную зависимость.

Спирмен мощнее Пирсона:

  • OUT_STAT/NEW10 - нелинейные данные,
  • R_B_COST/BOOK - выбросы,
  • log-ADD_FEE/OUT_STAT - выбросы,
  • BOOK/GRADUAT - выбросы.

3.3.1.2 Public:

Пирсон мощнее Спирмена:

  • R_B_COST/log-ADD_FEE - выбросы,
  • NEW10/log-ADD_FEE - выбросы усиливают линейную зависимость.

Спирмен мощнее Пирсона:

  • OUT_STAT/log-ADD_FEE - выбросы,
  • AVRCOMB/NEW10 - выбросы уменьшают коэффициент корреляции Пирсона.

3.3.2 Теоретическое объяснение (сравнение коэффициента корреляции Пирсона и Спирмена):

В тех ячейках, где совпадают или близки коэффициенты корреляции Пирсона и Спирмена почти нет выбросов и у scatterplot эллиптическая структура, если Пирсон значительно больше Спирмена - есть заметные выбросы. Спирмен может быть больше Пирсона, если выбросы имеются, но их положение не очень значимо относительно других наблюдений.

Значительные корреляции:

3.3.3 Private:

  • Выраженный блок GRADUAT, PH_D, SF_RATIO, AVRCOMB, NEW10, OUT_STAT.

3.3.4 Public:

  • Блок NEW10, AVRCOMB, PH_D, GRADUAT;
  • Пара NEW10, R_B_COST - нельзя сказать (по смыслу) о зависимости переменных. Могут зависеть от внешней переменной (скрытого фактора).

3.4 Проинтерпретируйте найденные корреляции - можно ли сказать, что является причиной, что следствием. Если есть какая-то другая причина, которая влияет одновременно на оба признака (скрытый фактор), то попробуйте убрать его влияние с помощью частных корреляций.

3.4.1 Private:

3.4.1.1 R_B_COST, SF_RATIO (hidden OUT_STAT)

Для R_B_COST и SF_RATIO выберем следующие скрытые факторы:

  • OUT_STAT - стоимость обучения связана с ценой проживания R_B_COST и качеством образования, что влияет на SF_RATIO.

Поставим уровень значимости alpha = 0.05

res<-pcor(col_split_mixed$private[,c("R_B_COST", "SF_RATIO", "OUT_STAT")])
res$estimate["R_B_COST", "SF_RATIO"]
## [1] -0.01433506
res$p.value["R_B_COST", "SF_RATIO"]
## [1] 0.9196585

estimate - получившаяся частная корреляция. Замечаем, что без OUT_STAT коэффициент корреляции Пирсона близок к 0, а p-value ~ 0.90.

3.4.1.2 AVRCOMB/BOOK (hidden OUT_STAT, NEW10)

  • Для AVRCOMB/BOOK можно выделить OUT_STAT по тем же причинам, что и в предыдущем примере и NEW10, так как они влияют на AVRCOMB
group_1<-c("AVRCOMB", "BOOK")
res<-pcor(na.omit(col_split_mixed$private[,c(group_1, "NEW10")]))
res
## $estimate
##           AVRCOMB        BOOK       NEW10
## AVRCOMB 1.0000000  0.15177454  0.95475775
## BOOK    0.1517745  1.00000000 -0.09376246
## NEW10   0.9547577 -0.09376246  1.00000000
## 
## $p.value
##              AVRCOMB      BOOK        NEW10
## AVRCOMB 0.000000e+00 0.4498169 1.122411e-14
## BOOK    4.498169e-01 0.0000000 6.418045e-01
## NEW10   1.122411e-14 0.6418045 0.000000e+00
## 
## $statistic
##            AVRCOMB       BOOK      NEW10
## AVRCOMB  0.0000000  0.7677672 16.0525714
## BOOK     0.7677672  0.0000000 -0.4708868
## NEW10   16.0525714 -0.4708868  0.0000000
## 
## $n
## [1] 28
## 
## $gp
## [1] 1
## 
## $method
## [1] "pearson"
res_2<-pcor(na.omit(col_split_mixed$private[,c(group_1, "OUT_STAT", "NEW10")]))
res_2
## $estimate
##            AVRCOMB        BOOK   OUT_STAT      NEW10
## AVRCOMB  1.0000000  0.23482510  0.2604398 0.81612411
## BOOK     0.2348251  1.00000000 -0.3815662 0.02045576
## OUT_STAT 0.2604398 -0.38156623  1.0000000 0.27217396
## NEW10    0.8161241  0.02045576  0.2721740 1.00000000
## 
## $p.value
##               AVRCOMB       BOOK   OUT_STAT        NEW10
## AVRCOMB  0.000000e+00 0.24820290 0.19879665 3.728780e-07
## BOOK     2.482029e-01 0.00000000 0.05442782 9.209917e-01
## OUT_STAT 1.987966e-01 0.05442782 0.00000000 1.785854e-01
## NEW10    3.728780e-07 0.92099170 0.17858544 0.000000e+00
## 
## $statistic
##           AVRCOMB       BOOK  OUT_STAT     NEW10
## AVRCOMB  0.000000  1.1834967  1.321494 6.9187347
## BOOK     1.183497  0.0000000 -2.022288 0.1002333
## OUT_STAT 1.321494 -2.0222884  0.000000 1.3856870
## NEW10    6.918735  0.1002333  1.385687 0.0000000
## 
## $n
## [1] 28
## 
## $gp
## [1] 2
## 
## $method
## [1] "pearson"

Получаем, что частная корреляция AVRCOMB/BOOK больше зависит от NEW10, чем от NEW10/OUT_STAT.

3.4.2 Public:

3.4.2.1 NEW10, R_B_COST (hidden OUT_STAT)

OUT_STAT - можно предположить, что цена обучения влияет на NEW10, R_B_COST.

group<-c("NEW10", "R_B_COST")
res<-pcor(na.omit(col_split_mixed$public[,c(group,"OUT_STAT")]), method="pearson")
res
## $estimate
##                NEW10  R_B_COST    OUT_STAT
## NEW10     1.00000000 0.5203072 -0.02803864
## R_B_COST  0.52030719 1.0000000  0.32506178
## OUT_STAT -0.02803864 0.3250618  1.00000000
## 
## $p.value
##                 NEW10     R_B_COST     OUT_STAT
## NEW10    0.000000e+00 9.240296e-09 0.7743545088
## R_B_COST 9.240296e-09 0.000000e+00 0.0006349232
## OUT_STAT 7.743545e-01 6.349232e-04 0.0000000000
## 
## $statistic
##               NEW10 R_B_COST   OUT_STAT
## NEW10     0.0000000 6.243199 -0.2874236
## R_B_COST  6.2431993 0.000000  3.5221714
## OUT_STAT -0.2874236 3.522171  0.0000000
## 
## $n
## [1] 108
## 
## $gp
## [1] 1
## 
## $method
## [1] "pearson"

Из estimate и p.value - отвергаем OUT_STAT как скрытый фактор.

3.4.2.2 GRADUAT, R_B_COST (hidden OUT_STAT, NEW10)

  • Зависимость между GRADUAT и R_B_COST. Скрытые факторы: OUT_STAT и NEW10.
group_2<-c("GRADUAT", "R_B_COST")
res_1<-pcor(na.omit(col_split_mixed$public[,c(group_2, "OUT_STAT", "NEW10")]), method = "spearman")
res_1
## $estimate
##            GRADUAT  R_B_COST   OUT_STAT      NEW10
## GRADUAT  1.0000000 0.1228687  0.4375959  0.4975827
## R_B_COST 0.1228687 1.0000000  0.2919737  0.2007210
## OUT_STAT 0.4375959 0.2919737  1.0000000 -0.2953932
## NEW10    0.4975827 0.2007210 -0.2953932  1.0000000
## 
## $p.value
##               GRADUAT    R_B_COST     OUT_STAT        NEW10
## GRADUAT  0.000000e+00 0.216289869 3.797030e-06 8.929736e-08
## R_B_COST 2.162899e-01 0.000000000 2.766667e-03 4.205514e-02
## OUT_STAT 3.797030e-06 0.002766667 0.000000e+00 2.451802e-03
## NEW10    8.929736e-08 0.042055136 2.451802e-03 0.000000e+00
## 
## $statistic
##           GRADUAT R_B_COST  OUT_STAT     NEW10
## GRADUAT  0.000000 1.244243  4.890929  5.764986
## R_B_COST 1.244243 0.000000  3.067983  2.059128
## OUT_STAT 4.890929 3.067983  0.000000 -3.107327
## NEW10    5.764986 2.059128 -3.107327  0.000000
## 
## $n
## [1] 105
## 
## $gp
## [1] 2
## 
## $method
## [1] "spearman"

Коэффициент корреляции Спирмена понизился, p-value ~ 0.22. Не отвергаем влияния совокупности скрытых факторов OUT_STAT и NEW10.

3.4.2.3 OUT_STAT, PH_D (hidden NEW10, SF_RATIO)

  • Коэффициент корреляции Спирмена OUT_STAT и PH_D равен 0.24, хотя визуально кажется, что они коррелируют сильнее. Скрытыми факторами могут быть NEW10, так как NEW10 отражает престижность вузов (NEW10 может зависеть от PH_D, поскольку лучшие студенты хотят найти лучших преподавателей), а OUT_STAT зависит от NEW10 и PH_D.
group_3<-c("PH_D", "OUT_STAT")
res_3<-pcor(na.omit(col_split_mixed$public[,c(group_3, "NEW10")]), method = "spearman")
res_3
## $estimate
##               PH_D    OUT_STAT       NEW10
## PH_D     1.0000000  0.16610803  0.50812933
## OUT_STAT 0.1661080  1.00000000 -0.04604289
## NEW10    0.5081293 -0.04604289  1.00000000
## 
## $p.value
##                  PH_D   OUT_STAT        NEW10
## PH_D     0.000000e+00 0.09520736 4.987348e-08
## OUT_STAT 9.520736e-02 0.00000000 6.458585e-01
## NEW10    4.987348e-08 0.64585848 0.000000e+00
## 
## $statistic
##              PH_D   OUT_STAT      NEW10
## PH_D     0.000000  1.6844818  5.8996943
## OUT_STAT 1.684482  0.0000000 -0.4609178
## NEW10    5.899694 -0.4609178  0.0000000
## 
## $n
## [1] 103
## 
## $gp
## [1] 1
## 
## $method
## [1] "spearman"

Не отвергаем влияние NEW10 на PH_D и OUT_STAT. Можно также рассмотреть частную корреляцию от SF_RATIO, так как это отражает (неявно) процент PH_D и стоимость OUT_STAT за счёт студентов и профессоров.

res_4<-pcor(na.omit(col_split_mixed$public[,c(group_3, "SF_RATIO")]), method = "spearman")
res_4
## $estimate
##                PH_D    OUT_STAT    SF_RATIO
## PH_D      1.0000000  0.20532911 -0.27979699
## OUT_STAT  0.2053291  1.00000000 -0.09968463
## SF_RATIO -0.2797970 -0.09968463  1.00000000
## 
## $p.value
##                 PH_D   OUT_STAT    SF_RATIO
## PH_D     0.000000000 0.02702694 0.002351145
## OUT_STAT 0.027026942 0.00000000 0.287028174
## SF_RATIO 0.002351145 0.28702817 0.000000000
## 
## $statistic
##               PH_D  OUT_STAT  SF_RATIO
## PH_D      0.000000  2.240044 -3.111698
## OUT_STAT  2.240044  0.000000 -1.069669
## SF_RATIO -3.111698 -1.069669  0.000000
## 
## $n
## [1] 117
## 
## $gp
## [1] 1
## 
## $method
## [1] "spearman"

Отвергаем влияние SF_RATIO как скрытого фактора. Не отвергаем, что PH_D скрытый фактор относительно SF_RATIO/OUT_STAT.

4 Регрессионный анализ

Задача: провести регрессию качества поступающих студентов (NEW10) на остальные признаки. С помощью пошаговой регрессии уменьшить число признаков

library(lm.beta)

4.1 Подготовка данных

4.1.1 Выбор признаков

Рассматриваемые наблюдения: 1. PPIND (для группировки) 2. ADD_FEE - дополнительные выплаты 3. BOOK - плата за бронирование 4. NEW10 (зависимая переменная, также NEW25) - студентов из топ 10% своей старшей школы 5. AVRCOMB - средний комбинированный балл 6. SF_RATIO - Student/Faculty ratio 7. PH_D - Pct. of faculty with Ph.D.’s 8. GRADUAT - Graduation rate 9. R_B_COST - Room and board costs 10. OUT_STAT - Out-of-state tuition

colNames<-c("...1", "ADD_FEE", "BOOK", "R_B_COST","OUT_STAT","NEW10","AVRCOMB","SF_RATIO","PH_D","GRADUAT","PPIND")
col_I_regr<-col_I_sn[c(colNames)]
col_I_regr<-col_I_regr|>mutate(ADD_FEE = log(ADD_FEE))|>rename(log_ADD_FEE=ADD_FEE)
col_I_regr_split<-split(col_I_regr, col_I_regr$PPIND)
# Автор: -Николай-
library(GGally)
library(ggplot2)
library(dplyr)
library(rlang)
## 
## Присоединяю пакет: 'rlang'
## Следующие объекты скрыты от 'package:purrr':
## 
##     %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
##     flatten_raw, invoke, splice
plot_numeric_pairs <- function(df, title = NULL, highlight = NULL) {
  # Выбираем имена всех numeric-столбцов
  numeric_cols <- df %>% dplyr::select(dplyr::where(is.numeric)) %>% names()
  
  # Функция для нижних панелей с возможностью выделения точек разными цветами
  custom_points <- function(data, mapping, ...) {
    # Отрисовываем базовые точки синим цветом
    p <- ggplot(data, mapping) +
      geom_point(alpha = 0.6, size = 2, color = "blue", ...)
    
    if (!is.null(highlight)) {
      # Если highlight не является списком, превращаем его в список с одним элементом и цветом "red"
      if (!is.list(highlight) || inherits(highlight, "data.frame")) {
        highlight <- list(red = highlight)
      }
      
      # Получаем имена переменных из mapping с помощью as_label()
      xvar <- as_label(mapping$x)
      yvar <- as_label(mapping$y)
      
      # Для каждого набора точек в списке highlight
      for (col in names(highlight)) {
        hl_df <- highlight[[col]]
        # Проверяем, что hl_df содержит нужные переменные
        if(all(c(xvar, yvar) %in% names(hl_df))) {
          # Находим совпадающие строки между данными панели и hl_df
          data_highlight <- semi_join(data, hl_df, by = c(xvar, yvar))
          # Накладываем выделенные точки нужного цвета
          p <- p + geom_point(data = data_highlight, mapping = mapping,
                              alpha = 0.6, size = 2, color = col, ...)
        }
      }
    }
    p
  }
  
  # Строим матрицу парных графиков
  ggpairs(df,
          columns      = numeric_cols,
          columnLabels = numeric_cols,
          title        = title,
          upper        = list(continuous = wrap("cor", size = 3)),
          lower        = list(continuous = custom_points),
          diag         = list(continuous = wrap("barDiag", bins = 10, fill = "lightblue", color = "black"))
  ) +
    theme_bw() +
    theme(
      strip.text.x  = element_text(size = 7, angle = 0, hjust = 0),
      strip.text.y  = element_text(size = 8, angle = 0, hjust = 0),
      axis.text = element_text(size = 6),
      plot.margin   = unit(c(1, 2, 1, 1), "cm")
    )
}

4.2 Анализ частных вузов

pr_data<-col_I_regr_split$private
print(paste0("Число наблюдений: ", dim(pr_data)[1]))
## [1] "Число наблюдений: 53"
colSums(is.na(pr_data))
##        ...1 log_ADD_FEE        BOOK    R_B_COST    OUT_STAT       NEW10 
##           0           7           0           0           0           1 
##     AVRCOMB    SF_RATIO        PH_D     GRADUAT       PPIND 
##          24           0           4           1           0

Много пропусков в AVRCOMB

plot_numeric_pairs(pr_data, "Частные вузы")

4.2.1 Выбросы “на глаз”

outlier_pr_book_eye<-pr_data[pr_data$BOOK>875,]
outlier_pr_out_eye<-pr_data[pr_data$OUT_STAT<5000,]
outlier_pr_book_eye
outlier_pr_out_eye
plot_numeric_pairs(pr_data, "Частные вузы", highlight = list(red=outlier_pr_book_eye, black=outlier_pr_out_eye))

4.2.2 Построение линейной регрессии

Здесь условия из .tsk файла

Стоит помнить, что у AVRCOMB много NA. В такой модели не будет выбросов “на глаз”

model_pr_new <- lm(NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT,
           data = pr_data)|>lm.beta()

4.2.3 Интерпретация результатов

summary(model_pr_new)
## 
## Call:
## lm(formula = NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + 
##     PH_D + SF_RATIO + GRADUAT + OUT_STAT, data = pr_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.4148 -3.8409  0.0271  2.9131 12.4176 
## 
## Coefficients:
##               Estimate Standardized Std. Error t value Pr(>|t|)  
## (Intercept) -78.090327           NA  51.790742  -1.508   0.1575  
## AVRCOMB       0.117625     0.544789   0.042417   2.773   0.0169 *
## BOOK         -0.006699    -0.029617   0.014003  -0.478   0.6410  
## log_ADD_FEE   3.200735     0.098598   2.219759   1.442   0.1749  
## R_B_COST     -0.003437    -0.117762   0.002499  -1.376   0.1941  
## PH_D         -0.492675    -0.126378   0.496416  -0.992   0.3406  
## SF_RATIO     -1.353492    -0.205756   0.594526  -2.277   0.0419 *
## GRADUAT       0.306142     0.178421   0.258574   1.184   0.2594  
## OUT_STAT      0.002354     0.257544   0.001391   1.693   0.1163  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.797 on 12 degrees of freedom
##   (32 пропущенных наблюдений удалены)
## Multiple R-squared:  0.9646, Adjusted R-squared:  0.9411 
## F-statistic: 40.92 on 8 and 12 DF,  p-value: 1.496e-07

4.2.4 Построение доверительных эллипсоидов

# Автор: -Вадим-
create_ellipse <- function(center, cov_matrix, level = 0.95, n = 100) {
  # Уровень доверия преобразуем в квантиль хи-квадрат
  chi_sq <- qchisq(level, df = 2)

  # Вычисляем собственные значения и собственные векторы
  eigen_res <- eigen(cov_matrix)

  # Углы для точек эллипса
  angles <- seq(0, 2 * pi, length.out = n)

  # Радиусы эллипса (основаны на собственных значениях)
  radii <- sqrt(chi_sq * eigen_res$values)
  # Создаем эллипс в стандартной ориентации
  ellipse_points <- cbind(radii[1] * cos(angles), radii[2] * sin(angles))

  # Поворачиваем эллипс в соответствии с собственными векторами
  ellipse_points <- t(eigen_res$vectors %*% t(ellipse_points))

  # Добавляем смещение
  ellipse_points <- sweep(ellipse_points, 2, center, "+")

  ellipse_points <- as.data.frame(ellipse_points)
  colnames(ellipse_points) <- c("x", "y") # Assign column names 

  return(ellipse_points)
}

plot_confidence_ellipse <- function(model, which.coeff = c(2, 4), level = 0.95, 
                               xlab = NULL, ylab = NULL, main = "Доверительный эллипс",
                               margin = 0.1) {
  
  # Проверка на то, что переданы именно 2 коэффициента
  if (length(which.coeff) != 2) {
    stop("which.coeff должен содержать ровно 2 индекса коэффициентов")
  }
  
  # Получение необходимых данных из модели
  model_summary <- summary(model)
  coef_estimates <- coef(model)
  сov_matrix <- model_summary$cov.unscaled
  
  # Проверка допустимости индексов коэффициентов
  if (any(which.coeff > length(coef_estimates))) {
    stop("Индекс коэффициента превышает их количество в модели")
  }
  
  # Создание эллипса
  ellipse_points <- create_ellipse(
    center = coef_estimates[which.coeff], 
    cov_matrix = сov_matrix[which.coeff, which.coeff], 
    level = level
  )
  
  # Определение границ графика с учетом отступов
  x_min <- min(ellipse_points$x)
  x_max <- max(ellipse_points$x)
  y_min <- min(ellipse_points$y)
  y_max <- max(ellipse_points$y)
  
  x_range <- c(x_min - margin * (x_max - x_min), x_max + margin * (x_max - x_min))
  y_range <- c(y_min - margin * (y_max - y_min), y_max + margin * (y_max - y_min))
  
  # Если не указаны названия осей, используем имена коэффициентов
  if (is.null(xlab)) {
    xlab <- names(coef_estimates)[which.coeff[1]]
  }
  if (is.null(ylab)) {
    ylab <- names(coef_estimates)[which.coeff[2]]
  }
  
  # Создаем data.frame с центральной точкой
  center_point <- data.frame(
    x = coef_estimates[which.coeff[1]],
    y = coef_estimates[which.coeff[2]]
  )
  
  # Создаем объект ggplot
  p <- ggplot2::ggplot() +
    # Добавляем эллипс
    ggplot2::geom_path(data = ellipse_points, ggplot2::aes(x = x, y = y), color = "blue") +
    # Добавляем точку с оценкой коэффициента
    ggplot2::geom_point(data = center_point, ggplot2::aes(x = x, y = y), color = "red", size = 3) +
    # Добавляем вертикальную и горизонтальную линии через центр
    ggplot2::geom_vline(xintercept = center_point$x, linetype = "dashed", color = "gray") +
    ggplot2::geom_hline(yintercept = center_point$y, linetype = "dashed", color = "gray") +
    # Настраиваем оси и заголовок
    ggplot2::xlim(x_range) +
    ggplot2::ylim(y_range) +
    ggplot2::labs(
      x = xlab,
      y = ylab,
      title = main
    ) +
    ggplot2::theme_light()
  
  return(p)
}
plot_confidence_ellipse(model_pr_new, c(3, 2), level=0.95)

На примере этой пары: AVRCOMB и SF_RATIO значимы

4.2.5 Влияние корреляции на качество оценки

На примере, PH_D и OUT_STAT

phcor<-cor(pr_data$PH_D, pr_data$OUT_STAT, use = "complete.obs")

y <- pr_data$NEW10
X <- pr_data[, c("PH_D", "OUT_STAT")]

model_cor <- lm(NEW10 ~ PH_D + OUT_STAT, data = pr_data)
sum_cor<-summary(model_cor)
sum_cor
## 
## Call:
## lm(formula = NEW10 ~ PH_D + OUT_STAT, data = pr_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.696 -13.720   0.979  15.266  34.815 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -1.009e+02  3.392e+01  -2.975  0.00470 **
## PH_D         1.551e+00  4.763e-01   3.257  0.00214 **
## OUT_STAT     1.198e-03  9.551e-04   1.255  0.21608   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 17.96 on 45 degrees of freedom
##   (5 пропущенных наблюдений удалены)
## Multiple R-squared:  0.4428, Adjusted R-squared:  0.4181 
## F-statistic: 17.88 on 2 and 45 DF,  p-value: 1.926e-06
# Оцениваем дисперсию ошибок
sigma_s_squared <- (sum_cor$sigma/sd(y, na.rm=TRUE))^2 #  оценка дисперсии ошибок
# Размер выборки
n <- nobs(model_cor)

R_XX<-cor(X, use="complete.obs")
print("Inverse R_XX")
## [1] "Inverse R_XX"
solve(R_XX)
##               PH_D  OUT_STAT
## PH_D      2.088331 -1.507579
## OUT_STAT -1.507579  2.088331
cov_beta<-sigma_s_squared/n*solve(R_XX)
print("Cov_beta:")
## [1] "Cov_beta:"
cov_beta
##                 PH_D    OUT_STAT
## PH_D      0.02360939 -0.01704377
## OUT_STAT -0.01704377  0.02360939
print("Determinant")
## [1] "Determinant"
detcoef<-1/(1-phcor^2)
detcoef
## [1] 2.088331
print("-rho/(1-rho^2)")
## [1] "-rho/(1-rho^2)"
-phcor*detcoef
## [1] -1.507579
phcor
## [1] 0.7219063

4.2.6 Множественная корреляция

# Автор: -Евгений-

# Вычисление множественной корреляции
calc_multicollinearity <- function(var, data) {
  X <- as.matrix(data[, setdiff(names(data), var)])  # Матрица остальных регрессоров
  y <- as.matrix(data[[var]])  # Текущий регрессор как зависимая переменная
  model <- lm(y ~ X) 
  sqrt(summary(model)$r.squared)  
}

calc_partial_correlation <- function(var, data) {
  # Остатки регрессии текущего регрессора на остальные регрессоры
  data<-data|>na.omit()
  model_X <- lm(data[[var]] ~ ., data = data[, !names(data) %in% c(var, "NEW10")])
  res_X <- residuals(model_X)
  
  # Остатки регрессии Y на остальные регрессоры
  model_Y <- lm(NEW10 ~ ., data = data[, !names(data) %in% var])
  res_Y <- residuals(model_Y)
  
  cor(res_X, res_Y, use = "complete.obs")  # Корреляция остатков
}

regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
  Регрессор = regressors,
  Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
  Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)

print(redundancy_table)
##               Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE                0.6803712          0.3842868
## BOOK               BOOK                0.4955889         -0.1368050
## R_B_COST       R_B_COST                0.8079355         -0.3690491
## OUT_STAT       OUT_STAT                0.9472357          0.4390291
## AVRCOMB         AVRCOMB                0.9764603          0.6249356
## SF_RATIO       SF_RATIO                0.8649231         -0.5492085
## PH_D               PH_D                0.9121799         -0.2754193
## GRADUAT         GRADUAT                0.9401219          0.3234121

Сильно положительные множественные корреляции вместе с частными корреляциями, следует сравнить по отдельности

regr_high <- c("R_B_COST", "OUT_STAT", "AVRCOMB", "SF_RATIO", "PH_D", "GRADUAT")
regr_fullhigh<-c(regr_high, "NEW10")
redund_high <- data.frame(
  Регрессор = regr_high,
  Множественная_корреляция = sapply(regr_high, calc_multicollinearity, pr_data[, regr_fullhigh]),
  Частная_корреляция = sapply(regr_high, calc_partial_correlation, pr_data[,regr_fullhigh])
)

print(redund_high)
##          Регрессор Множественная_корреляция Частная_корреляция
## R_B_COST  R_B_COST                0.6621454         -0.1634685
## OUT_STAT  OUT_STAT                0.9214250          0.3142483
## AVRCOMB    AVRCOMB                0.9687653          0.5706967
## SF_RATIO  SF_RATIO                0.8369497         -0.3982266
## PH_D          PH_D                0.8830917         -0.1158525
## GRADUAT    GRADUAT                0.9332980          0.3403546

Видно, что PH_D имеет высокую множественную корреляцию и низкую частную корреляцию, поэтому выбросим её

model_pr_noPhD<-lm(NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + SF_RATIO + GRADUAT + OUT_STAT, data=pr_data)|>lm.beta()
summary(model_pr_noPhD)
## 
## Call:
## lm(formula = NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + 
##     SF_RATIO + GRADUAT + OUT_STAT, data = pr_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.9908 -2.5824  0.6409  3.1699 11.8038 
## 
## Coefficients:
##               Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept) -1.156e+02           NA  2.833e+01  -4.080 0.000872 ***
## AVRCOMB      1.291e-01    6.067e-01  2.746e-02   4.703 0.000240 ***
## BOOK        -4.336e-03   -1.885e-02  1.215e-02  -0.357 0.725883    
## log_ADD_FEE  2.568e+00    7.431e-02  2.010e+00   1.277 0.219712    
## R_B_COST    -3.258e-03   -1.042e-01  2.189e-03  -1.488 0.156128    
## SF_RATIO    -1.300e+00   -1.876e-01  5.055e-01  -2.572 0.020479 *  
## GRADUAT      1.590e-01    8.799e-02  1.657e-01   0.960 0.351347    
## OUT_STAT     1.837e-03    1.893e-01  1.079e-03   1.704 0.107803    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.753 on 16 degrees of freedom
##   (29 пропущенных наблюдений удалены)
## Multiple R-squared:  0.9616, Adjusted R-squared:  0.9448 
## F-statistic: 57.28 on 7 and 16 DF,  p-value: 3.83e-10
AIC(model_pr_noPhD)
## [1] 168.0601
regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10", "PH_D"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
  Регрессор = regressors,
  Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
  Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)

print(redundancy_table)
##               Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE                0.5973852         0.30419786
## BOOK               BOOK                0.3840774        -0.08885404
## R_B_COST       R_B_COST                0.7551732        -0.34870910
## OUT_STAT       OUT_STAT                0.9141040         0.39183509
## AVRCOMB         AVRCOMB                0.9692789         0.76171345
## SF_RATIO       SF_RATIO                0.8252002        -0.54081295
## GRADUAT         GRADUAT                0.8544490         0.23336928

Попробуем убрать BOOK (из-за низкой частной корреляции) и GRADUAT

model_pr1<-lm(NEW10 ~ AVRCOMB + log_ADD_FEE + R_B_COST + SF_RATIO + OUT_STAT, data=pr_data)|>lm.beta()
summary(model_pr1)
## 
## Call:
## lm(formula = NEW10 ~ AVRCOMB + log_ADD_FEE + R_B_COST + SF_RATIO + 
##     OUT_STAT, data = pr_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.4607  -1.9482  -0.7267   2.3280  13.4017 
## 
## Coefficients:
##               Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept) -1.185e+02           NA  2.745e+01  -4.315 0.000417 ***
## AVRCOMB      1.372e-01    6.446e-01  2.323e-02   5.907 1.36e-05 ***
## log_ADD_FEE  2.961e+00    8.567e-02  1.880e+00   1.575 0.132718    
## R_B_COST    -3.886e-03   -1.243e-01  2.040e-03  -1.905 0.072923 .  
## SF_RATIO    -1.336e+00   -1.927e-01  4.914e-01  -2.718 0.014095 *  
## OUT_STAT     2.164e-03    2.229e-01  9.870e-04   2.192 0.041756 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.583 on 18 degrees of freedom
##   (29 пропущенных наблюдений удалены)
## Multiple R-squared:  0.959,  Adjusted R-squared:  0.9476 
## F-statistic: 84.15 on 5 and 18 DF,  p-value: 7.676e-12
AIC(model_pr1)
## [1] 165.6631
regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10", "PH_D", "BOOK", "GRADUAT"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
  Регрессор = regressors,
  Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
  Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)

print(redundancy_table)
##               Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE                0.5684235          0.3479827
## R_B_COST       R_B_COST                0.7447506         -0.4095647
## OUT_STAT       OUT_STAT                0.9088596          0.4590401
## AVRCOMB         AVRCOMB                0.9668878          0.8122173
## SF_RATIO       SF_RATIO                0.8236772         -0.5394659
AIC(model_pr_new)
## [1] 148.3376
AIC(model_pr_noPhD)
## [1] 168.0601
AIC(model_pr1)
## [1] 165.6631

4.2.7 OLS. Пошаговый AIC.

library(olsrr)
## 
## Присоединяю пакет: 'olsrr'
## Следующий объект скрыт от 'package:MASS':
## 
##     cement
## Следующий объект скрыт от 'package:datasets':
## 
##     rivers
backward_aic<-ols_step_backward_aic(model_pr_new, progress=TRUE, details=TRUE)
## Backward Elimination Method 
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. AVRCOMB 
## 2. BOOK 
## 3. log_ADD_FEE 
## 4. R_B_COST 
## 5. PH_D 
## 6. SF_RATIO 
## 7. GRADUAT 
## 8. OUT_STAT 
## 
## 
## Step     => 0 
## Model    => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT 
## AIC      => 148.3376 
## 
## Initiating stepwise selection... 
## 
##                   Table: Removing Existing Variables                    
## -----------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC       R2       Adj. R2 
## -----------------------------------------------------------------------
## BOOK            1    146.734    156.135    97.586    0.96397    0.94457 
## PH_D            1    147.994    157.395    97.513    0.96174    0.94114 
## GRADUAT         1    148.658    158.058    97.491    0.96051    0.93925 
## R_B_COST        1    149.412    158.813    97.480    0.95907    0.93703 
## log_ADD_FEE     1    149.693    159.094    97.480    0.95852    0.93618 
## OUT_STAT        1    150.834    160.235    97.501    0.95620    0.93262 
## SF_RATIO        1    153.877    163.277    97.737    0.94937    0.92211 
## AVRCOMB         1    156.737    166.137    98.197    0.94199    0.91075 
## -----------------------------------------------------------------------
## 
## Step     => 1 
## Removed  => BOOK 
## Model    => NEW10 ~ AVRCOMB + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT 
## AIC      => 146.7343 
## 
##                   Table: Removing Existing Variables                    
## -----------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC       R2       Adj. R2 
## -----------------------------------------------------------------------
## PH_D            1    147.170    155.526    94.423    0.95954    0.94220 
## GRADUAT         1    147.861    156.217    94.521    0.95818    0.94026 
## log_ADD_FEE     1    147.962    156.319    94.536    0.95798    0.93997 
## R_B_COST        1    148.265    156.622    94.583    0.95737    0.93910 
## OUT_STAT        1    150.933    159.289    95.079    0.95160    0.93085 
## SF_RATIO        1    153.017    161.373    95.578    0.94655    0.92364 
## AVRCOMB         1    155.056    163.412    96.162    0.94110    0.91585 
## -----------------------------------------------------------------------
## 
## 
## No more variables to be removed.
## 
## Variables Removed: 
## 
## => BOOK
forward_aic<-ols_step_forward_aic(model_pr_new, progress=TRUE, details=TRUE)
## Forward Selection Method 
## ------------------------
## 
## Candidate Terms: 
## 
## 1. AVRCOMB 
## 2. BOOK 
## 3. log_ADD_FEE 
## 4. R_B_COST 
## 5. PH_D 
## 6. SF_RATIO 
## 7. GRADUAT 
## 8. OUT_STAT 
## 
## 
## Step     => 0 
## Model    => NEW10 ~ 1 
## AIC      => 202.5251 
## 
## Initiating stepwise selection... 
## 
##                        Table: Adding New Variables                        
## -------------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC        R2       Adj. R2  
## -------------------------------------------------------------------------
## AVRCOMB         1    147.617    150.751     87.731    0.93346     0.92995 
## OUT_STAT        1    171.142    174.276    107.789    0.79601     0.78527 
## GRADUAT         1    176.803    179.937    112.954    0.73289     0.71883 
## PH_D            1    186.429    189.563    121.963    0.57756     0.55533 
## SF_RATIO        1    189.574    192.707    124.955    0.50933     0.48350 
## log_ADD_FEE     1    202.476    205.609    137.417    0.09297     0.04523 
## R_B_COST        1    203.410    206.543    138.328    0.05173     0.00182 
## BOOK            1    204.456    207.590    139.350    0.00327    -0.04919 
## -------------------------------------------------------------------------
## 
## Step     => 1 
## Added    => AVRCOMB 
## Model    => NEW10 ~ AVRCOMB 
## AIC      => 147.617 
## 
##                       Table: Adding New Variables                       
## -----------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC       R2       Adj. R2 
## -----------------------------------------------------------------------
## SF_RATIO        1    145.049    149.227    86.340    0.94647    0.94052 
## OUT_STAT        1    147.493    151.671    88.069    0.93986    0.93318 
## BOOK            1    147.975    152.153    88.412    0.93846    0.93162 
## log_ADD_FEE     1    148.581    152.759    88.846    0.93666    0.92962 
## PH_D            1    149.095    153.274    89.215    0.93509    0.92788 
## GRADUAT         1    149.166    153.344    89.265    0.93487    0.92764 
## R_B_COST        1    149.275    153.453    89.344    0.93453    0.92726 
## -----------------------------------------------------------------------
## 
## Step     => 2 
## Added    => SF_RATIO 
## Model    => NEW10 ~ AVRCOMB + SF_RATIO 
## AIC      => 145.0489 
## 
##                       Table: Adding New Variables                       
## -----------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC       R2       Adj. R2 
## -----------------------------------------------------------------------
## BOOK            1    145.266    150.488    87.603    0.95082    0.94215 
## OUT_STAT        1    145.397    150.620    87.679    0.95052    0.94178 
## GRADUAT         1    146.151    151.373    88.121    0.94871    0.93966 
## log_ADD_FEE     1    146.479    151.701    88.314    0.94790    0.93871 
## R_B_COST        1    146.891    152.113    88.557    0.94687    0.93749 
## PH_D            1    146.968    152.190    88.603    0.94667    0.93726 
## -----------------------------------------------------------------------
## 
## 
## No more variables to be added.
## 
## Variables Selected: 
## 
## => AVRCOMB 
## => SF_RATIO

Получили, что forward_aic имеет лучшее значение, чем backward_aic, то есть, является более хорошей линейной моделью. К тому же, у forward_aic меньше признаков.

4.2.7.1 Forward

forward_aic
## 
## 
##                              Stepwise Summary                              
## -------------------------------------------------------------------------
## Step    Variable        AIC        SBC       SBIC        R2       Adj. R2 
## -------------------------------------------------------------------------
##  0      Base Model    202.525    204.614    139.293    0.00000    0.00000 
##  1      AVRCOMB       147.617    150.751     87.731    0.93346    0.92995 
##  2      SF_RATIO      145.049    149.227     86.340    0.94647    0.94052 
## -------------------------------------------------------------------------
## 
## Final Model Output 
## ------------------
## 
##                          Model Summary                          
## ---------------------------------------------------------------
## R                       0.973       RMSE                 6.323 
## R-Squared               0.946       MSE                 39.975 
## Adj. R-Squared          0.941       Coef. Var           11.813 
## Pred R-Squared          0.925       AIC                145.049 
## MAE                     4.914       SBC                149.227 
## ---------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
##  AIC: Akaike Information Criteria 
##  SBC: Schwarz Bayesian Criteria 
## 
##                                 ANOVA                                  
## ----------------------------------------------------------------------
##                  Sum of                                               
##                 Squares        DF    Mean Square       F         Sig. 
## ----------------------------------------------------------------------
## Regression    14841.759         2       7420.879    159.117    0.0000 
## Residual        839.479        18         46.638                      
## Total         15681.238        20                                     
## ----------------------------------------------------------------------
## 
##                                      Parameter Estimates                                       
## ----------------------------------------------------------------------------------------------
##       model        Beta    Std. Error    Std. Beta      t        Sig        lower       upper 
## ----------------------------------------------------------------------------------------------
## (Intercept)    -150.963        21.538                 -7.009    0.000    -196.212    -105.714 
##     AVRCOMB       0.188         0.015        0.869    12.124    0.000       0.155       0.220 
##    SF_RATIO      -0.986         0.471       -0.150    -2.091    0.051      -1.977       0.004 
## ----------------------------------------------------------------------------------------------

4.2.7.2 Backward

backward_aic
## 
## 
##                              Stepwise Summary                              
## -------------------------------------------------------------------------
## Step    Variable        AIC        SBC       SBIC        R2       Adj. R2 
## -------------------------------------------------------------------------
##  0      Full Model    148.338    158.783    101.117    0.96464    0.94107 
##  1      BOOK          146.734    156.135     96.228    0.96397    0.94457 
## -------------------------------------------------------------------------
## 
## Final Model Output 
## ------------------
## 
##                          Model Summary                          
## ---------------------------------------------------------------
## R                       0.982       RMSE                 5.187 
## R-Squared               0.964       MSE                 26.905 
## Adj. R-Squared          0.945       Coef. Var           11.404 
## Pred R-Squared          0.901       AIC                146.734 
## MAE                     4.081       SBC                156.135 
## ---------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
##  AIC: Akaike Information Criteria 
##  SBC: Schwarz Bayesian Criteria 
## 
##                                 ANOVA                                 
## ---------------------------------------------------------------------
##                  Sum of                                              
##                 Squares        DF    Mean Square      F         Sig. 
## ---------------------------------------------------------------------
## Regression    15116.223         7       2159.460    49.685    0.0000 
## Residual        565.015        13         43.463                     
## Total         15681.238        20                                    
## ---------------------------------------------------------------------
## 
##                                     Parameter Estimates                                     
## -------------------------------------------------------------------------------------------
##       model       Beta    Std. Error    Std. Beta      t        Sig        lower     upper 
## -------------------------------------------------------------------------------------------
## (Intercept)    -71.623        48.490                 -1.477    0.163    -176.379    33.133 
##     AVRCOMB      0.110         0.038        0.511     2.873    0.013       0.027     0.193 
## log_ADD_FEE      3.162         2.151        0.097     1.470    0.165      -1.486     7.810 
##    R_B_COST     -0.004         0.002       -0.126    -1.543    0.147      -0.009     0.001 
##        PH_D     -0.573         0.453       -0.147    -1.264    0.228      -1.552     0.406 
##    SF_RATIO     -1.413         0.564       -0.215    -2.507    0.026      -2.631    -0.195 
##     GRADUAT      0.344         0.238        0.201     1.445    0.172      -0.171     0.860 
##    OUT_STAT      0.003         0.001        0.286     2.113    0.055       0.000     0.005 
## -------------------------------------------------------------------------------------------

Здесь можно сравнить по R^2, например.

model_pr_new<-lm(NEW10 ~ AVRCOMB + SF_RATIO, data=pr_data)|>lm.beta()
summary(model_pr_new)
## 
## Call:
## lm(formula = NEW10 ~ AVRCOMB + SF_RATIO, data = pr_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.0779  -6.4604   0.1359   4.5401  17.2921 
## 
## Coefficients:
##               Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept) -133.57926           NA   20.45055  -6.532 7.66e-07 ***
## AVRCOMB        0.17283      0.84064    0.01434  12.055 6.52e-12 ***
## SF_RATIO      -1.18802     -0.17520    0.47288  -2.512   0.0188 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.591 on 25 degrees of freedom
##   (25 пропущенных наблюдений удалены)
## Multiple R-squared:  0.9309, Adjusted R-squared:  0.9254 
## F-statistic: 168.4 on 2 and 25 DF,  p-value: 3.104e-15
AIC(model_pr_new)
## [1] 197.7974

4.2.8 Построение и анализ графиков

shapiro.test(residuals(model_pr_new))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_pr_new)
## W = 0.9525, p-value = 0.2288
plot(model_pr_new, 1:5)

У 15-го сильное плечо (ещё видно, как линия идёт к ней)

4.2.9 Поиск outliers по Residuals vs. Deleted Residuals

plot(residuals(model_pr_new), rstudent(model_pr_new))

Все значения не так сильно отличаются, так что посмотрим на outliers по Куку и Махаланобису

4.2.10 Поиск outliers по Cook и Mahalanobis

4.2.10.1 Cook

# Рассчитаем Cook's distance для модели с логарифмами
cooks_vals <- cooks.distance(model_pr_new)
cook_df <- data.frame(
  obs  = seq_along(cooks_vals),
  cook = cooks_vals
) %>%
  arrange(desc(cook))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(cook_df))


plot(
  x, cook_df$cook,
  type  = "h",
  main  = "Cook's Distance",
  xlab  = "Номер наблюдения (из исходного df)",
  ylab  = "Cook's distance",
  xaxt  = "n"                # убираем автоматическую ось X
)
axis(side = 1, at = x, labels = cook_df$obs, las = 2, cex.axis = 0.8)
abline(h = 1, col = "red", lty = 2)

4.2.10.2 Mahalanobis

# 7.2. Расстояние Махаланобиса (Mahalanobis Distance)
# Расстояние Махаланобиса в пространстве предикторов (регрессоров)
mahalanobis_distance <- mahalanobis(model.matrix(model_pr_new)[,-1],
                                     center = colMeans(model.matrix(model_pr_new)[,-1]),
                                     cov = cov(model.matrix(model_pr_new)[,-1]))
mah_df <- data.frame(
  obs  = seq_along(mahalanobis_distance),
  mah = mahalanobis_distance
) %>%
  arrange(desc(mah))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(mah_df))
alpha<-0.05

# Scatter plot для расстояния Махаланобиса
plot(
  x, mah_df$mah,
  type  = "h",
  main  = "Mahalanobis's Distance",
  xlab  = "Номер наблюдения (из исходного df)",
  ylab  = "Mahalanobis's distance",
  xaxt  = "n"                
)
axis(side = 1, at = x, labels = mah_df$obs, las = 2, cex.axis = 0.8)
abline(h = qchisq(p = 1 - alpha, df = nrow(mah_df)), col = "red", lty = 2)

Не нашёл outliers (в случае AVRCOMB)

4.2.11 Итог: прогнозирование

set.seed(456) # Другое зерно для генерации новых значений
max_AVRCOMB<-max(pr_data$AVRCOMB, na.rm=TRUE)
max_SF_RATIO<-max(pr_data$SF_RATIO, na.rm=TRUE)
new_data <- data.frame(
  AVRCOMB = max_AVRCOMB + 1:8,
  SF_RATIO = max_SF_RATIO + 1:8
)

# Прогнозируем значение NEW10
predicted_values <- predict(model_pr_new, newdata = new_data)
cat("Прогнозируемые значения NEW10:\n")
## Прогнозируемые значения NEW10:
print(predicted_values)
##        1        2        3        4        5        6        7        8 
## 84.73871 83.72351 82.70832 81.69312 80.67793 79.66273 78.64754 77.63234
# 6. Доверительные и предсказательные интервалы

# Доверительный интервал (confidence interval) для среднего значения NEW10
confidence_intervals <- predict(model_pr_new, newdata = new_data, interval = "confidence", level = 0.95)
print("Доверительные интервалы (95%):\n")
## [1] "Доверительные интервалы (95%):\n"
print(confidence_intervals)
##        fit      lwr      upr
## 1 84.73871 68.04653 101.4309
## 2 83.72351 66.08992 101.3571
## 3 82.70832 64.12789 101.2887
## 4 81.69312 62.16121 101.2250
## 5 80.67793 60.19053 101.1653
## 6 79.66273 58.21640 101.1091
## 7 78.64754 56.23924 101.0558
## 8 77.63234 54.25944 101.0052
# Предсказательный интервал (prediction interval) для отдельного наблюдения NEW10
prediction_intervals <- predict(model_pr_new, newdata = new_data, interval = "prediction", level = 0.95)
print("Предсказательные интервалы (95%):\n")
## [1] "Предсказательные интервалы (95%):\n"
print(prediction_intervals)
##        fit      lwr      upr
## 1 84.73871 61.86841 107.6090
## 2 83.72351 60.15733 107.2897
## 3 82.70832 58.42552 106.9911
## 4 81.69312 56.67478 106.7115
## 5 80.67793 54.90672 106.4491
## 6 79.66273 53.12282 106.2026
## 7 78.64754 51.32441 105.9707
## 8 77.63234 49.51270 105.7520
# Создаем data.frame для прогнозов и интервалов
predictions_df <- data.frame(
  AVRCOMB = new_data$AVRCOMB,
  SF_RATIO = new_data$SF_RATIO,
  Predicted = predicted_values,
  Lower_Conf = confidence_intervals[, "lwr"],
  Upper_Conf = confidence_intervals[, "upr"],
  Lower_Pred = prediction_intervals[, "lwr"],
  Upper_Pred = prediction_intervals[, "upr"]
)
predictions_df
# Визуализация (только для демонстрации, нужен более сложный график для нескольких предикторов)
ggplot(predictions_df, aes(x = 1:nrow(predictions_df), y = Predicted)) + # Условная ось x, т.к. несколько предикторов
  geom_point() +
  geom_errorbar(aes(ymin = Lower_Conf, ymax = Upper_Conf), color = "blue", width = 0.2) + # Доверительные интервалы
  geom_errorbar(aes(ymin = Lower_Pred, ymax = Upper_Pred), color = "red", width = 0.2) + # Предсказательные интервалы
  labs(title = "Прогнозы NEW10 с доверительными и предсказательными интервалами",
       x = "Прогноз (номер)",
       y = "NEW10") +
  scale_x_continuous(breaks = 1:nrow(predictions_df)) + # Подписи для оси x
  theme_bw()

  # annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.9, label = "Синие: Доверительные интервалы", color = "blue", hjust = 0) +
  # annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.8, label = "Красные: Предсказательные интервалы", color = "red", hjust = 0)

4.3 Анализ общественных вузов

pub_data<-col_I_regr_split$public
print(paste0("Число наблюдений: ", dim(pub_data)[1]))
## [1] "Число наблюдений: 123"
colSums(is.na(pub_data))
##        ...1 log_ADD_FEE        BOOK    R_B_COST    OUT_STAT       NEW10 
##           0          33           2           0           1          15 
##     AVRCOMB    SF_RATIO        PH_D     GRADUAT       PPIND 
##          43           0           5           4           0

Много пропусков в AVRCOMB

plot_numeric_pairs(pub_data, "Общественные вузы")

4.3.1 Выбросы “на глаз”

outlier_pub_out_eye<-pub_data[pub_data$OUT_STAT>14001,]
outlier_pub_out_eye
plot_numeric_pairs(pub_data, "Общественные вузы", highlight = list(red=outlier_pub_out_eye))

4.3.2 Построение линейной регрессии

Здесь условия из .tsk файла

Стоит помнить, что у AVRCOMB много NA. В такой модели не будет выбросов “на глаз”

model_pub_new <- lm(NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT,
           data = pub_data)|>lm.beta()

4.3.3 Интерпретация результатов

summary(model_pub_new)
## 
## Call:
## lm(formula = NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + 
##     PH_D + SF_RATIO + GRADUAT + OUT_STAT, data = pub_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.462  -7.954  -1.811   5.800  33.454 
## 
## Coefficients:
##               Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept) -1.957e+02           NA  5.205e+01  -3.760 0.000531 ***
## AVRCOMB      1.124e-01    3.828e-01  4.176e-02   2.691 0.010254 *  
## BOOK         3.210e-02    1.521e-01  2.409e-02   1.333 0.190033    
## log_ADD_FEE  4.645e+00    2.355e-01  2.416e+00   1.922 0.061536 .  
## R_B_COST     4.859e-03    1.906e-01  3.400e-03   1.429 0.160534    
## PH_D         4.937e-01    1.328e-01  4.894e-01   1.009 0.318999    
## SF_RATIO     7.921e-02    1.507e-02  6.456e-01   0.123 0.902959    
## GRADUAT      1.657e-01    1.096e-01  2.419e-01   0.685 0.497348    
## OUT_STAT    -8.391e-04   -9.592e-02  1.337e-03  -0.628 0.533709    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.48 on 41 degrees of freedom
##   (73 пропущенных наблюдений удалены)
## Multiple R-squared:  0.5774, Adjusted R-squared:  0.4949 
## F-statistic: 7.002 on 8 and 41 DF,  p-value: 8.678e-06

4.3.4 Построение доверительных эллипсоидов

plot_confidence_ellipse(model_pub_new, c(3, 2), level=0.95)

На примере этой пары: AVRCOMB и SF_RATIO значимы

4.3.5 Множественная корреляция

regressors <- setdiff(names(pub_data), c("...1", "PPIND", "NEW10"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
  Регрессор = regressors,
  Множественная_корреляция = sapply(regressors, calc_multicollinearity, pub_data[, regressors_full]),
  Частная_корреляция = sapply(regressors, calc_partial_correlation, pub_data[,regressors_full])
)

print(redundancy_table)
##               Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE                0.6083961          0.2875337
## BOOK               BOOK                0.4915332          0.2037467
## R_B_COST       R_B_COST                0.6691166          0.2178369
## OUT_STAT       OUT_STAT                0.7502290         -0.0975559
## AVRCOMB         AVRCOMB                0.7530387          0.3874765
## SF_RATIO       SF_RATIO                0.5634165          0.0191561
## PH_D               PH_D                0.6475179          0.1556248
## GRADUAT         GRADUAT                0.7761031          0.1063354

Можем убрать OUT_STAT и GRADUAT

AIC(model_pub_new)
## [1] 425.9309
model_pub1 <- lm(NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO,
           data = pub_data)|>lm.beta()
summary(model_pub1)
## 
## Call:
## lm(formula = NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + 
##     PH_D + SF_RATIO, data = pub_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.863  -7.446  -2.303   7.075  32.472 
## 
## Coefficients:
##               Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept) -2.035e+02           NA  4.107e+01  -4.955 1.07e-05 ***
## AVRCOMB      1.151e-01    4.081e-01  3.088e-02   3.727 0.000539 ***
## BOOK         3.398e-02    1.619e-01  2.134e-02   1.592 0.118345    
## log_ADD_FEE  5.201e+00    2.632e-01  2.298e+00   2.263 0.028499 *  
## R_B_COST     3.886e-03    1.546e-01  2.844e-03   1.367 0.178554    
## PH_D         5.734e-01    1.550e-01  4.594e-01   1.248 0.218431    
## SF_RATIO     1.030e-01    1.985e-02  5.417e-01   0.190 0.849989    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.14 on 45 degrees of freedom
##   (71 пропущенное наблюдение удалено)
## Multiple R-squared:  0.5601, Adjusted R-squared:  0.5014 
## F-statistic: 9.548 on 6 and 45 DF,  p-value: 9.133e-07
AIC(model_pub1)
## [1] 438.6636

AIC увеличился, а значит предпочтительнее выбирать model_pub_new.

regr_high <- setdiff(names(pub_data), c("...1", "PPIND", "NEW10", "OUT_STAT", "GRADUAT"))
regr_fullhigh<-c(regr_high, "NEW10")
redund_high <- data.frame(
  Регрессор = regr_high,
  Множественная_корреляция = sapply(regr_high, calc_multicollinearity, pub_data[, regr_fullhigh]),
  Частная_корреляция = sapply(regr_high, calc_partial_correlation, pub_data[,regr_fullhigh])
)

print(redund_high)
##               Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE                0.5924855         0.31966785
## BOOK               BOOK                0.3246989         0.23093274
## R_B_COST       R_B_COST                0.5163902         0.19961448
## AVRCOMB         AVRCOMB                0.6139954         0.48569728
## SF_RATIO       SF_RATIO                0.3206426         0.02834558
## PH_D               PH_D                0.6222610         0.18292398

4.3.6 OLS. Пошаговый AIC.

library(olsrr)
backward_aic_pub<-ols_step_backward_aic(model_pub_new, progress=TRUE, details=TRUE)
## Backward Elimination Method 
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. AVRCOMB 
## 2. BOOK 
## 3. log_ADD_FEE 
## 4. R_B_COST 
## 5. PH_D 
## 6. SF_RATIO 
## 7. GRADUAT 
## 8. OUT_STAT 
## 
## 
## Step     => 0 
## Model    => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT 
## AIC      => 425.9309 
## 
## Initiating stepwise selection... 
## 
##                    Table: Removing Existing Variables                    
## ------------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC        R2       Adj. R2 
## ------------------------------------------------------------------------
## SF_RATIO        1    423.949    441.157    285.464    0.57725    0.50679 
## OUT_STAT        1    424.409    441.617    285.755    0.57334    0.50223 
## GRADUAT         1    424.499    441.708    285.813    0.57257    0.50133 
## PH_D            1    425.157    442.365    286.230    0.56691    0.49473 
## BOOK            1    426.051    443.259    286.802    0.55910    0.48562 
## R_B_COST        1    426.362    443.570    287.002    0.55635    0.48241 
## log_ADD_FEE     1    428.246    445.454    288.223    0.53932    0.46253 
## AVRCOMB         1    432.065    449.273    290.751    0.50275    0.41987 
## ------------------------------------------------------------------------
## 
## Step     => 1 
## Removed  => SF_RATIO 
## Model    => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + GRADUAT + OUT_STAT 
## AIC      => 423.9492 
## 
##                    Table: Removing Existing Variables                    
## ------------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC        R2       Adj. R2 
## ------------------------------------------------------------------------
## GRADUAT         1    422.552    437.849    283.437    0.57212    0.51241 
## OUT_STAT        1    422.682    437.978    283.525    0.57101    0.51115 
## PH_D            1    423.231    438.527    283.901    0.56627    0.50575 
## BOOK            1    424.053    439.349    284.466    0.55908    0.49755 
## R_B_COST        1    424.548    439.844    284.807    0.55470    0.49256 
## log_ADD_FEE     1    426.379    441.676    286.079    0.53808    0.47363 
## AVRCOMB         1    430.065    445.361    288.678    0.50275    0.43336 
## ------------------------------------------------------------------------
## 
## Step     => 2 
## Removed  => GRADUAT 
## Model    => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D + OUT_STAT 
## AIC      => 422.5524 
## 
##                    Table: Removing Existing Variables                    
## ------------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC        R2       Adj. R2 
## ------------------------------------------------------------------------
## OUT_STAT        1    420.913    434.297    281.292    0.56902    0.52004 
## PH_D            1    421.997    435.381    282.087    0.55957    0.50952 
## BOOK            1    422.446    435.830    282.417    0.55560    0.50510 
## R_B_COST        1    423.127    436.511    282.920    0.54951    0.49831 
## log_ADD_FEE     1    425.924    439.308    284.998    0.52358    0.46945 
## AVRCOMB         1    435.531    448.915    292.326    0.42266    0.35706 
## ------------------------------------------------------------------------
## 
## Step     => 3 
## Removed  => OUT_STAT 
## Model    => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST + PH_D 
## AIC      => 420.9132 
## 
##                    Table: Removing Existing Variables                    
## ------------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC        R2       Adj. R2 
## ------------------------------------------------------------------------
## PH_D            1    420.513    431.985    280.150    0.55501    0.51545 
## R_B_COST        1    421.160    432.632    280.657    0.54921    0.50914 
## BOOK            1    421.831    433.303    281.185    0.54312    0.50250 
## log_ADD_FEE     1    424.596    436.068    283.366    0.51715    0.47423 
## AVRCOMB         1    433.721    445.193    290.696    0.42046    0.36895 
## ------------------------------------------------------------------------
## 
## Step     => 4 
## Removed  => PH_D 
## Model    => NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST 
## AIC      => 420.5128 
## 
##                    Table: Removing Existing Variables                    
## ------------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC        R2       Adj. R2 
## ------------------------------------------------------------------------
## R_B_COST        1    421.792    431.352    280.561    0.52484    0.49386 
## BOOK            1    422.106    431.666    280.822    0.52185    0.49067 
## log_ADD_FEE     1    425.733    435.293    283.858    0.48588    0.45235 
## AVRCOMB         1    436.611    446.171    293.093    0.36093    0.31925 
## ------------------------------------------------------------------------
## 
## 
## No more variables to be removed.
## 
## Variables Removed: 
## 
## => SF_RATIO 
## => GRADUAT 
## => OUT_STAT 
## => PH_D
forward_aic_pub<-ols_step_forward_aic(model_pub_new, progress=TRUE, details=TRUE)
## Forward Selection Method 
## ------------------------
## 
## Candidate Terms: 
## 
## 1. AVRCOMB 
## 2. BOOK 
## 3. log_ADD_FEE 
## 4. R_B_COST 
## 5. PH_D 
## 6. SF_RATIO 
## 7. GRADUAT 
## 8. OUT_STAT 
## 
## 
## Step     => 0 
## Model    => NEW10 ~ 1 
## AIC      => 452.9976 
## 
## Initiating stepwise selection... 
## 
##                        Table: Adding New Variables                        
## -------------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC        R2       Adj. R2  
## -------------------------------------------------------------------------
## AVRCOMB         1    432.770    438.506    290.015    0.35889     0.34553 
## PH_D            1    440.325    446.061    297.005    0.25432     0.23878 
## log_ADD_FEE     1    440.669    446.405    297.324    0.24917     0.23353 
## GRADUAT         1    440.775    446.511    297.423    0.24757     0.23189 
## R_B_COST        1    445.964    451.700    302.247    0.16529     0.14790 
## BOOK            1    450.887    456.623    306.843    0.07893     0.05974 
## OUT_STAT        1    454.461    460.197    310.192    0.01067    -0.00994 
## SF_RATIO        1    454.634    460.370    310.354    0.00724    -0.01344 
## -------------------------------------------------------------------------
## 
## Step     => 1 
## Added    => AVRCOMB 
## Model    => NEW10 ~ AVRCOMB 
## AIC      => 432.7703 
## 
##                       Table: Adding New Variables                        
## ------------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC        R2       Adj. R2 
## ------------------------------------------------------------------------
## log_ADD_FEE     1    423.754    431.402    281.872    0.48566    0.46377 
## R_B_COST        1    426.389    434.038    284.194    0.45782    0.43475 
## PH_D            1    427.487    435.135    285.163    0.44579    0.42221 
## BOOK            1    431.895    439.543    289.066    0.39471    0.36895 
## GRADUAT         1    433.164    440.812    290.193    0.37916    0.35274 
## OUT_STAT        1    434.743    442.391    291.599    0.35923    0.33197 
## SF_RATIO        1    434.747    442.395    291.602    0.35919    0.33192 
## ------------------------------------------------------------------------
## 
## Step     => 2 
## Added    => log_ADD_FEE 
## Model    => NEW10 ~ AVRCOMB + log_ADD_FEE 
## AIC      => 423.7541 
## 
##                      Table: Adding New Variables                       
## ----------------------------------------------------------------------
## Predictor    DF      AIC        SBC       SBIC        R2       Adj. R2 
## ----------------------------------------------------------------------
## BOOK          1    421.792    431.352    280.561    0.52484    0.49386 
## R_B_COST      1    422.106    431.666    280.822    0.52185    0.49067 
## PH_D          1    422.181    431.741    280.885    0.52113    0.48990 
## OUT_STAT      1    425.569    435.129    283.721    0.48756    0.45414 
## GRADUAT       1    425.579    435.140    283.729    0.48745    0.45403 
## SF_RATIO      1    425.746    435.306    283.869    0.48574    0.45221 
## ----------------------------------------------------------------------
## 
## Step     => 3 
## Added    => BOOK 
## Model    => NEW10 ~ AVRCOMB + log_ADD_FEE + BOOK 
## AIC      => 421.7919 
## 
##                      Table: Adding New Variables                       
## ----------------------------------------------------------------------
## Predictor    DF      AIC        SBC       SBIC        R2       Adj. R2 
## ----------------------------------------------------------------------
## R_B_COST      1    420.513    431.985    280.150    0.55501    0.51545 
## PH_D          1    421.160    432.632    280.657    0.54921    0.50914 
## GRADUAT       1    422.984    434.456    282.092    0.53246    0.49091 
## OUT_STAT      1    423.748    435.220    282.695    0.52526    0.48306 
## SF_RATIO      1    423.783    435.255    282.723    0.52493    0.48270 
## ----------------------------------------------------------------------
## 
## Step     => 4 
## Added    => R_B_COST 
## Model    => NEW10 ~ AVRCOMB + log_ADD_FEE + BOOK + R_B_COST 
## AIC      => 420.5128 
## 
##                      Table: Adding New Variables                       
## ----------------------------------------------------------------------
## Predictor    DF      AIC        SBC       SBIC        R2       Adj. R2 
## ----------------------------------------------------------------------
## PH_D          1    420.913    434.297    281.292    0.56902    0.52004 
## OUT_STAT      1    421.997    435.381    282.087    0.55957    0.50952 
## GRADUAT       1    422.228    435.613    282.257    0.55753    0.50725 
## SF_RATIO      1    422.498    435.882    282.455    0.55514    0.50459 
## ----------------------------------------------------------------------
## 
## 
## No more variables to be added.
## 
## Variables Selected: 
## 
## => AVRCOMB 
## => log_ADD_FEE 
## => BOOK 
## => R_B_COST

Получили совпадение результатов пошагового forward и пошагового backward.

4.3.6.1 Forward

forward_aic_pub
## 
## 
##                               Stepwise Summary                              
## --------------------------------------------------------------------------
## Step    Variable         AIC        SBC       SBIC        R2       Adj. R2 
## --------------------------------------------------------------------------
##  0      Base Model     452.998    456.822    309.665    0.00000    0.00000 
##  1      AVRCOMB        432.770    438.506    290.015    0.35889    0.34553 
##  2      log_ADD_FEE    423.754    431.402    281.872    0.48566    0.46377 
##  3      BOOK           421.792    431.352    280.561    0.52484    0.49386 
##  4      R_B_COST       420.513    431.985    280.150    0.55501    0.51545 
## --------------------------------------------------------------------------
## 
## Final Model Output 
## ------------------
## 
##                          Model Summary                           
## ----------------------------------------------------------------
## R                        0.745       RMSE                14.385 
## R-Squared                0.555       MSE                206.930 
## Adj. R-Squared           0.515       Coef. Var           42.498 
## Pred R-Squared           0.424       AIC                420.513 
## MAE                     11.265       SBC                431.985 
## ----------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
##  AIC: Akaike Information Criteria 
##  SBC: Schwarz Bayesian Criteria 
## 
##                                 ANOVA                                 
## ---------------------------------------------------------------------
##                  Sum of                                              
##                 Squares        DF    Mean Square      F         Sig. 
## ---------------------------------------------------------------------
## Regression    12904.403         4       3226.101    14.031    0.0000 
## Residual      10346.477        45        229.922                     
## Total         23250.880        49                                    
## ---------------------------------------------------------------------
## 
##                                      Parameter Estimates                                       
## ----------------------------------------------------------------------------------------------
##       model        Beta    Std. Error    Std. Beta      t        Sig        lower       upper 
## ----------------------------------------------------------------------------------------------
## (Intercept)    -184.284        31.816                 -5.792    0.000    -248.364    -120.203 
##     AVRCOMB       0.136         0.031        0.462     4.430    0.000       0.074       0.197 
## log_ADD_FEE       5.845         2.211        0.296     2.644    0.011       1.392      10.297 
##        BOOK       0.039         0.021        0.185     1.831    0.074      -0.004       0.082 
##    R_B_COST       0.005         0.003        0.192     1.746    0.088      -0.001       0.011 
## ----------------------------------------------------------------------------------------------

4.3.6.2 Backward

backward_aic_pub
## 
## 
##                              Stepwise Summary                              
## -------------------------------------------------------------------------
## Step    Variable        AIC        SBC       SBIC        R2       Adj. R2 
## -------------------------------------------------------------------------
##  0      Full Model    425.931    445.051    287.892    0.57740    0.49494 
##  1      SF_RATIO      423.949    441.157    285.030    0.57725    0.50679 
##  2      GRADUAT       422.552    437.849    282.885    0.57212    0.51241 
##  3      OUT_STAT      420.913    434.297    280.618    0.56902    0.52004 
##  4      PH_D          420.513    431.985    279.705    0.55501    0.51545 
## -------------------------------------------------------------------------
## 
## Final Model Output 
## ------------------
## 
##                          Model Summary                           
## ----------------------------------------------------------------
## R                        0.745       RMSE                14.385 
## R-Squared                0.555       MSE                206.930 
## Adj. R-Squared           0.515       Coef. Var           42.498 
## Pred R-Squared           0.424       AIC                420.513 
## MAE                     11.265       SBC                431.985 
## ----------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
##  AIC: Akaike Information Criteria 
##  SBC: Schwarz Bayesian Criteria 
## 
##                                 ANOVA                                 
## ---------------------------------------------------------------------
##                  Sum of                                              
##                 Squares        DF    Mean Square      F         Sig. 
## ---------------------------------------------------------------------
## Regression    12904.403         4       3226.101    14.031    0.0000 
## Residual      10346.477        45        229.922                     
## Total         23250.880        49                                    
## ---------------------------------------------------------------------
## 
##                                      Parameter Estimates                                       
## ----------------------------------------------------------------------------------------------
##       model        Beta    Std. Error    Std. Beta      t        Sig        lower       upper 
## ----------------------------------------------------------------------------------------------
## (Intercept)    -184.284        31.816                 -5.792    0.000    -248.364    -120.203 
##     AVRCOMB       0.136         0.031        0.462     4.430    0.000       0.074       0.197 
##        BOOK       0.039         0.021        0.185     1.831    0.074      -0.004       0.082 
## log_ADD_FEE       5.845         2.211        0.296     2.644    0.011       1.392      10.297 
##    R_B_COST       0.005         0.003        0.192     1.746    0.088      -0.001       0.011 
## ----------------------------------------------------------------------------------------------
model_pub_new<-lm(NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST, data=pub_data)|>lm.beta()
summary(model_pub_new)
## 
## Call:
## lm(formula = NEW10 ~ AVRCOMB + BOOK + log_ADD_FEE + R_B_COST, 
##     data = pub_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.246  -9.548  -1.513   9.069  29.958 
## 
## Coefficients:
##               Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept) -1.895e+02           NA  2.957e+01  -6.410 5.02e-08 ***
## AVRCOMB      1.261e-01    4.023e-01  2.927e-02   4.308 7.70e-05 ***
## BOOK         4.740e-02    2.089e-01  2.025e-02   2.341   0.0233 *  
## log_ADD_FEE  6.998e+00    3.365e-01  2.135e+00   3.278   0.0019 ** 
## R_B_COST     5.708e-03    2.171e-01  2.666e-03   2.141   0.0371 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.23 on 50 degrees of freedom
##   (68 пропущенных наблюдений удалены)
## Multiple R-squared:  0.6167, Adjusted R-squared:  0.5861 
## F-statistic: 20.11 on 4 and 50 DF,  p-value: 6.357e-10
AIC(model_pub_new)
## [1] 462.4016

4.3.7 Построение и анализ графиков

shapiro.test(residuals(model_pub_new))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_pub_new)
## W = 0.97633, p-value = 0.3475
plot(model_pub_new)

На первом графике Residuals vs. Fitted заметен сильный изгиб (матожидание ошибки ненулевое)

4.3.8 Поиск outliers по Residuals vs. Deleted Residuals

plot(residuals(model_pub_new), rstudent(model_pub_new))

Все значения не так сильно отличаются, так что посмотрим на outliers по Куку и Махаланобису

4.3.9 Поиск outliers по Cook и Mahalanobis

4.3.9.1 Cook

# Рассчитаем Cook's distance для модели с логарифмами
cooks_vals <- cooks.distance(model_pub_new)
cook_df <- data.frame(
  obs  = seq_along(cooks_vals),
  cook = cooks_vals
) %>%
  arrange(desc(cook))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(cook_df))


plot(
  x, cook_df$cook,
  type  = "h",
  main  = "Cook's Distance",
  xlab  = "Номер наблюдения (из исходного df)",
  ylab  = "Cook's distance",
  xaxt  = "n"                # убираем автоматическую ось X
)
axis(side = 1, at = x, labels = cook_df$obs, las = 2, cex.axis = 0.8)
abline(h = 1, col = "red", lty = 2)

4.3.9.2 Mahalanobis

# 7.2. Расстояние Махаланобиса (Mahalanobis Distance)
# Расстояние Махаланобиса в пространстве предикторов (регрессоров)
mahalanobis_distance <- mahalanobis(model.matrix(model_pub_new)[,-1],
                                     center = colMeans(model.matrix(model_pub_new)[,-1]),
                                     cov = cov(model.matrix(model_pub_new)[,-1]))

mah_df <- data.frame(
  obs  = seq_along(mahalanobis_distance),
  mah = mahalanobis_distance
) %>%
  arrange(desc(mah))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(mah_df))
alpha<-0.05

# Scatter plot для расстояния Махаланобиса
plot(
  x, mah_df$mah,
  type  = "h",
  main  = "Mahalanobis's Distance",
  xlab  = "Номер наблюдения (из исходного df)",
  ylab  = "Mahalanobis's distance",
  xaxt  = "n"                
)
axis(side = 1, at = x, labels = mah_df$obs, las = 2, cex.axis = 0.8)
abline(h = qchisq(p = 1 - alpha, df = nrow(mah_df)), col = "red", lty = 2)

Не нашёл outliers (в случае AVRCOMB)

4.3.10 Итог: прогнозирование

set.seed(456) # Другое зерно для генерации новых значений
max_AVRCOMB<-max(pub_data$AVRCOMB, na.rm=TRUE)
max_BOOK<-max(pub_data$BOOK, na.rm=TRUE)
max_log_ADD_FEE<-max(pub_data$log_ADD_FEE, na.rm=TRUE)
max_R_B_COST<-max(pub_data$R_B_COST, na.rm=TRUE)
new_data <- data.frame(
  AVRCOMB = max_AVRCOMB - seq(50, 120, by = 10),
  BOOK = max_BOOK - seq(10, 80, by = 10),
  log_ADD_FEE = max_log_ADD_FEE - seq(0.01, 0.08, by=0.01),
  R_B_COST = max_R_B_COST - seq(10, 80, by=10)
)
# Прогнозируем значение NEW10
predicted_values <- predict(model_pub_new, newdata = new_data)
cat("Прогнозируемые значения NEW10:\n")
## Прогнозируемые значения NEW10:
print(predicted_values)
##        1        2        3        4        5        6        7        8 
## 96.97503 95.11299 93.25095 91.38891 89.52687 87.66483 85.80280 83.94076
# 6. Доверительные и предсказательные интервалы

# Доверительный интервал (confidence interval) для среднего значения NEW10
confidence_intervals <- predict(model_pub_new, newdata = new_data, interval = "confidence", level = 0.95)
print("Доверительные интервалы (95%):\n")
## [1] "Доверительные интервалы (95%):\n"
print(confidence_intervals)
##        fit      lwr       upr
## 1 96.97503 81.95568 111.99437
## 2 95.11299 80.48533 109.74065
## 3 93.25095 78.99581 107.50609
## 4 91.38891 77.48558 105.29224
## 5 89.52687 75.95304 103.10071
## 6 87.66483 74.39651 100.93315
## 7 85.80280 72.81432  98.79127
## 8 83.94076 71.20476  96.67675
# Предсказательный интервал (prediction interval) для отдельного наблюдения NEW10
prediction_intervals <- predict(model_pub_new, newdata = new_data, interval = "prediction", level = 0.95)
print("Предсказательные интервалы (95%):\n")
## [1] "Предсказательные интервалы (95%):\n"
print(prediction_intervals)
##        fit      lwr      upr
## 1 96.97503 62.89617 131.0539
## 2 95.11299 61.20493 129.0210
## 3 93.25095 59.50192 127.0000
## 4 91.38891 57.78697 124.9909
## 5 89.52687 56.05992 122.9938
## 6 87.66483 54.32063 121.0090
## 7 85.80280 52.56895 119.0366
## 8 83.94076 50.80477 117.0767
predictions_df <- data.frame(
  AVRCOMB = new_data$AVRCOMB,
  BOOK = new_data$BOOK,
  log_ADD_FEE = new_data$log_ADD_FEE,
  R_B_COST = new_data$R_B_COST,
  Predicted = predicted_values,
  Lower_Conf = confidence_intervals[, "lwr"],
  Upper_Conf = confidence_intervals[, "upr"],
  Lower_Pred = prediction_intervals[, "lwr"],
  Upper_Pred = prediction_intervals[, "upr"]
)
predictions_df
ggplot(predictions_df, aes(x = 1:nrow(predictions_df), y = Predicted)) + 
  geom_point() +
  geom_errorbar(aes(ymin = Lower_Conf, ymax = Upper_Conf), color = "blue", width = 0.2) + # Доверительные интервалы
  geom_errorbar(aes(ymin = Lower_Pred, ymax = Upper_Pred), color = "red", width = 0.2) + # Предсказательные интервалы
  labs(title = "Прогнозы NEW10 с доверительными и предсказательными интервалами",
       x = "Прогноз (номер)",
       y = "NEW10") +
  scale_x_continuous(breaks = 1:nrow(predictions_df)) + # Подписи для оси x
  theme_bw()

4.4 Исправления.

Важные пункты для исправления:

  • Использование AVRCOMB (особенно в случае частных вузов) не интересно. Поэтому уберём её из рассматриваемых признаков.

  • Доверительный эллипсоид некорректно строится, поэтому нужно использовать другую функцию построения доверительного эллипсоида.

  • Рассмотреть прологарифмированный NEW10 для общественных вузов

4.5 Подготовка данных

4.5.1 Выбор признаков

Рассматриваемые наблюдения: 1. PPIND (для группировки) 2. ADD_FEE - дополнительные выплаты 3. BOOK - плата за бронирование 4. NEW10 (зависимая переменная, также NEW25) - студентов из топ 10% своей старшей школы 5. AVRCOMB - средний комбинированный балл (лишний признак) 6. SF_RATIO - Student/Faculty ratio 7. PH_D - Pct. of faculty with Ph.D.’s 8. GRADUAT - Graduation rate 9. R_B_COST - Room and board costs 10. OUT_STAT - Out-of-state tuition

colNames<-c("...1", "ADD_FEE", "BOOK", "R_B_COST","OUT_STAT","NEW10", "SF_RATIO","PH_D","GRADUAT","PPIND")
col_I_regr<-col_I_sn[c(colNames)]
col_I_regr<-col_I_regr|>mutate(ADD_FEE = log(ADD_FEE))|>rename(log_ADD_FEE=ADD_FEE)
col_I_regr_split<-split(col_I_regr, col_I_regr$PPIND)
# Автор: -Николай-
library(GGally)
library(ggplot2)
library(dplyr)
library(rlang)

plot_numeric_pairs <- function(df, title = NULL, highlight = NULL) {
  # Выбираем имена всех numeric-столбцов
  numeric_cols <- df %>% dplyr::select(dplyr::where(is.numeric)) %>% names()
  
  # Функция для нижних панелей с возможностью выделения точек разными цветами
  custom_points <- function(data, mapping, ...) {
    # Отрисовываем базовые точки синим цветом
    p <- ggplot(data, mapping) +
      geom_point(alpha = 0.6, size = 2, color = "blue", ...)
    
    if (!is.null(highlight)) {
      # Если highlight не является списком, превращаем его в список с одним элементом и цветом "red"
      if (!is.list(highlight) || inherits(highlight, "data.frame")) {
        highlight <- list(red = highlight)
      }
      
      # Получаем имена переменных из mapping с помощью as_label()
      xvar <- as_label(mapping$x)
      yvar <- as_label(mapping$y)
      
      # Для каждого набора точек в списке highlight
      for (col in names(highlight)) {
        hl_df <- highlight[[col]]
        # Проверяем, что hl_df содержит нужные переменные
        if(all(c(xvar, yvar) %in% names(hl_df))) {
          # Находим совпадающие строки между данными панели и hl_df
          data_highlight <- semi_join(data, hl_df, by = c(xvar, yvar))
          # Накладываем выделенные точки нужного цвета
          p <- p + geom_point(data = data_highlight, mapping = mapping,
                              alpha = 0.6, size = 2, color = col, ...)
        }
      }
    }
    p
  }
  
  # Строим матрицу парных графиков
  ggpairs(df,
          columns      = numeric_cols,
          columnLabels = numeric_cols,
          title        = title,
          upper        = list(continuous = wrap("cor", size = 3)),
          lower        = list(continuous = custom_points),
          diag         = list(continuous = wrap("barDiag", bins = 10, fill = "lightblue", color = "black"))
  ) +
    theme_bw() +
    theme(
      strip.text.x  = element_text(size = 7, angle = 0, hjust = 0),
      strip.text.y  = element_text(size = 8, angle = 0, hjust = 0),
      axis.text = element_text(size = 6),
      plot.margin   = unit(c(1, 2, 1, 1), "cm")
    )
}

4.6 Общий парный график

ggpairs(col_I_regr[,-which(names(col_I_regr)=="...1")], aes(alpha = 0.5, color=PPIND),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Комментарий: замечаем, что для публичных вузов NEW10 имеет выраженный правый хвост. Поэтому стоит использовать log_NEW10 для публичных вузов.

При постановке вопроса о построении линейной модели с категориальным признаком стоит обратить внимание на то, близки ли корреляции разных групп. Если они отличаются, то лучше рассматривать каждую категорию по отдельности. Если близки, то вместе. Это связано с тем, что разные корреляции по разному влияют на общую модель, будут возникать неоднородности.

4.7 Анализ частных вузов

pr_data<-col_I_regr_split$private
print(paste0("Число наблюдений: ", dim(pr_data)[1]))
## [1] "Число наблюдений: 53"
colSums(is.na(pr_data))
##        ...1 log_ADD_FEE        BOOK    R_B_COST    OUT_STAT       NEW10 
##           0           7           0           0           0           1 
##    SF_RATIO        PH_D     GRADUAT       PPIND 
##           0           4           1           0

4.7.1 График

plot_numeric_pairs(pr_data, "Частные вузы")

4.7.2 Выбросы “на глаз”

outlier_pr_book_eye<-pr_data[pr_data$BOOK>875,]
outlier_pr_out_eye<-pr_data[pr_data$OUT_STAT<5000,]
outlier_pr_book_eye
outlier_pr_out_eye
plot_numeric_pairs(pr_data, "Частные вузы", list(red = outlier_pr_book_eye, green = outlier_pr_out_eye))

4.7.3 Линейная модель

library(lm.beta)
model_pr_full <- lm(NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT, 
                    data = pr_data) |> lm.beta()
summary(model_pr_full)
## 
## Call:
## lm(formula = NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + 
##     GRADUAT + OUT_STAT, data = pr_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -25.1636  -6.8538   0.0723   7.2414  25.5012 
## 
## Coefficients:
##              Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept) -9.289705           NA  38.856630  -0.239  0.81252    
## BOOK         0.002300     0.015211   0.013944   0.165  0.86999    
## log_ADD_FEE  1.877765     0.063416   2.978846   0.630  0.53280    
## R_B_COST    -0.002897    -0.120242   0.002971  -0.975  0.33658    
## PH_D        -0.219550    -0.069638   0.508017  -0.432  0.66843    
## SF_RATIO    -1.772400    -0.324548   0.604579  -2.932  0.00608 ** 
## GRADUAT      1.069124     0.642349   0.215538   4.960 2.07e-05 ***
## OUT_STAT     0.001308     0.192473   0.001098   1.192  0.24180    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.4 on 33 degrees of freedom
##   (12 пропущенных наблюдений удалены)
## Multiple R-squared:   0.73,  Adjusted R-squared:  0.6727 
## F-statistic: 12.75 on 7 and 33 DF,  p-value: 8.508e-08

4.7.4 Доверительные эллипсы

plot_confidence_ellipse(model_pr_full, c(2, 3))

plot_confidence_ellipse(model_pr_full, c(4, 5))

plot_confidence_ellipse(model_pr_full, c(6, 7))

plot_confidence_ellipse(model_pr_full, c(8, 2))

# Автор: -Евгений-

# Вычисление множественной корреляции
calc_multicollinearity <- function(var, data) {
  X <- as.matrix(data[, setdiff(names(data), var)])  # Матрица остальных регрессоров
  y <- as.matrix(data[[var]])  # Текущий регрессор как зависимая переменная
  model <- lm(y ~ X) 
  sqrt(summary(model)$r.squared)  
}

calc_partial_correlation <- function(var, data) {
  # Остатки регрессии текущего регрессора на остальные регрессоры
  data<-data|>na.omit()
  model_X <- lm(data[[var]] ~ ., data = data[, !names(data) %in% c(var, "NEW10")])
  res_X <- residuals(model_X)
  
  # Остатки регрессии Y на остальные регрессоры
  model_Y <- lm(NEW10 ~ ., data = data[, !names(data) %in% var])
  res_Y <- residuals(model_Y)
  
  cor(res_X, res_Y, use = "complete.obs")  # Корреляция остатков
}

regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
  Регрессор = regressors,
  Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
  Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)

print(redundancy_table)
##               Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE                0.4485260         0.10907796
## BOOK               BOOK                0.1965671         0.02870252
## R_B_COST       R_B_COST                0.6906122        -0.16735732
## OUT_STAT       OUT_STAT                0.8361740         0.20315778
## SF_RATIO       SF_RATIO                0.6858132        -0.45456011
## PH_D               PH_D                0.8286465        -0.07501937
## GRADUAT         GRADUAT                0.8488216         0.65354784

Комментарий: Стоит исключить PH_D из-за высокой множественной и низкой частной корреляций (также можно исключить BOOK из-за низкой частной корреляции).

AIC(model_pr_full)
## [1] 338.2841
summary(model_pr_full)
## 
## Call:
## lm(formula = NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + 
##     GRADUAT + OUT_STAT, data = pr_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -25.1636  -6.8538   0.0723   7.2414  25.5012 
## 
## Coefficients:
##              Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept) -9.289705           NA  38.856630  -0.239  0.81252    
## BOOK         0.002300     0.015211   0.013944   0.165  0.86999    
## log_ADD_FEE  1.877765     0.063416   2.978846   0.630  0.53280    
## R_B_COST    -0.002897    -0.120242   0.002971  -0.975  0.33658    
## PH_D        -0.219550    -0.069638   0.508017  -0.432  0.66843    
## SF_RATIO    -1.772400    -0.324548   0.604579  -2.932  0.00608 ** 
## GRADUAT      1.069124     0.642349   0.215538   4.960 2.07e-05 ***
## OUT_STAT     0.001308     0.192473   0.001098   1.192  0.24180    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.4 on 33 degrees of freedom
##   (12 пропущенных наблюдений удалены)
## Multiple R-squared:   0.73,  Adjusted R-squared:  0.6727 
## F-statistic: 12.75 on 7 and 33 DF,  p-value: 8.508e-08
model_pr_nono<-lm(NEW10 ~ log_ADD_FEE + R_B_COST + SF_RATIO + GRADUAT + OUT_STAT, data=pr_data)|>lm.beta()
summary(model_pr_nono)
## 
## Call:
## lm(formula = NEW10 ~ log_ADD_FEE + R_B_COST + SF_RATIO + GRADUAT + 
##     OUT_STAT, data = pr_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.1896  -7.8512   0.4492   6.5775  29.4475 
## 
## Coefficients:
##               Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept) -1.254e+01           NA  2.304e+01  -0.544   0.5895    
## log_ADD_FEE  1.313e+00    4.221e-02  2.810e+00   0.467   0.6430    
## R_B_COST    -3.726e-03   -1.462e-01  2.743e-03  -1.358   0.1823    
## SF_RATIO    -1.872e+00   -3.275e-01  5.701e-01  -3.284   0.0022 ** 
## GRADUAT      1.019e+00    5.877e-01  1.880e-01   5.420 3.56e-06 ***
## OUT_STAT     1.228e-03    1.710e-01  9.368e-04   1.311   0.1976    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.49 on 38 degrees of freedom
##   (9 пропущенных наблюдений удалены)
## Multiple R-squared:  0.726,  Adjusted R-squared:  0.6899 
## F-statistic: 20.13 on 5 and 38 DF,  p-value: 9.083e-10
AIC(model_pr_nono)
## [1] 361.386

Комментарий: замечаем, что после того, как убрали BOOK и PH_D увеличился Adj.R^2, немного уменьшился R^2, уменьшился p-value, но увеличился AIC.

regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10", "PH_D", "BOOK"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
  Регрессор = regressors,
  Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
  Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)

print(redundancy_table)
##               Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE                0.3480743         0.07558012
## R_B_COST       R_B_COST                0.6376041        -0.21521296
## OUT_STAT       OUT_STAT                0.7708944         0.20809285
## SF_RATIO       SF_RATIO                0.6597111        -0.47020518
## GRADUAT         GRADUAT                0.8087224         0.66027976

Попробуем убрать log_ADD_FEE

model_pr1<-lm(NEW10 ~ GRADUAT + R_B_COST + SF_RATIO + OUT_STAT, data=pr_data)|>lm.beta()
summary(model_pr1)
## 
## Call:
## lm(formula = NEW10 ~ GRADUAT + R_B_COST + SF_RATIO + OUT_STAT, 
##     data = pr_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -27.041  -9.947  -1.776   8.355  49.615 
## 
## Coefficients:
##               Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept) 14.1924930           NA 21.7789414   0.652  0.51786    
## GRADUAT      0.9698293    0.6097925  0.1999565   4.850 1.45e-05 ***
## R_B_COST    -0.0038847   -0.1556376  0.0029142  -1.333  0.18908    
## SF_RATIO    -1.6934700   -0.3080606  0.6201252  -2.731  0.00893 ** 
## OUT_STAT     0.0003467    0.0543364  0.0010089   0.344  0.73268    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.36 on 46 degrees of freedom
##   (2 пропущенных наблюдений удалены)
## Multiple R-squared:  0.6291, Adjusted R-squared:  0.5968 
## F-statistic:  19.5 on 4 and 46 DF,  p-value: 1.918e-09
AIC(model_pr1)
## [1] 430.1013

Комментарий: AIC увеличилась, а R^2 сильно уменьшились, p-value увеличился.

regressors <- setdiff(names(pr_data), c("...1", "PPIND", "NEW10", "PH_D", "BOOK", "log_ADD_FEE"))
regressors_full<-c(regressors, "NEW10")
redundancy_table <- data.frame(
  Регрессор = regressors,
  Множественная_корреляция = sapply(regressors, calc_multicollinearity, pr_data[, regressors_full]),
  Частная_корреляция = sapply(regressors, calc_partial_correlation, pr_data[,regressors_full])
)

print(redundancy_table)
##          Регрессор Множественная_корреляция Частная_корреляция
## R_B_COST  R_B_COST                0.6561089        -0.19285565
## OUT_STAT  OUT_STAT                0.8235937         0.05060282
## SF_RATIO  SF_RATIO                0.6743579        -0.37350259
## GRADUAT    GRADUAT                0.8139325         0.58168914
AIC(model_pr_full)
## [1] 338.2841
AIC(model_pr_nono)
## [1] 361.386
AIC(model_pr1)
## [1] 430.1013

4.7.5 OLS. Пошаговый AIC.

4.7.5.1 Backward

library(olsrr)
backward_aic<-ols_step_backward_aic(model_pr_full, progress=TRUE, details=TRUE)
## Backward Elimination Method 
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. BOOK 
## 2. log_ADD_FEE 
## 3. R_B_COST 
## 4. PH_D 
## 5. SF_RATIO 
## 6. GRADUAT 
## 7. OUT_STAT 
## 
## 
## Step     => 0 
## Model    => NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT 
## AIC      => 338.2841 
## 
## Initiating stepwise selection... 
## 
##                    Table: Removing Existing Variables                    
## ------------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC        R2       Adj. R2 
## ------------------------------------------------------------------------
## BOOK            1    336.318    350.026    223.228    0.72977    0.68209 
## PH_D            1    336.516    350.224    223.348    0.72847    0.68055 
## log_ADD_FEE     1    336.775    350.483    223.505    0.72674    0.67852 
## R_B_COST        1    337.449    351.157    223.916    0.72221    0.67319 
## OUT_STAT        1    338.012    351.721    224.262    0.71837    0.66867 
## SF_RATIO        1    345.774    359.483    229.220    0.65968    0.59962 
## GRADUAT         1    359.125    372.833    238.570    0.52868    0.44551 
## ------------------------------------------------------------------------
## 
## Step     => 1 
## Removed  => BOOK 
## Model    => NEW10 ~ log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT 
## AIC      => 336.3179 
## 
##                    Table: Removing Existing Variables                    
## ------------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC        R2       Adj. R2 
## ------------------------------------------------------------------------
## PH_D            1    334.535    346.530    220.890    0.72834    0.68953 
## log_ADD_FEE     1    334.821    346.816    221.080    0.72644    0.68736 
## R_B_COST        1    335.500    347.495    221.536    0.72187    0.68213 
## OUT_STAT        1    336.045    348.040    221.902    0.71815    0.67789 
## SF_RATIO        1    343.786    355.781    227.258    0.65958    0.61094 
## GRADUAT         1    357.284    369.279    237.270    0.52684    0.45925 
## ------------------------------------------------------------------------
## 
## Step     => 2 
## Removed  => PH_D 
## Model    => NEW10 ~ log_ADD_FEE + R_B_COST + SF_RATIO + GRADUAT + OUT_STAT 
## AIC      => 334.5351 
## 
##                    Table: Removing Existing Variables                    
## ------------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC        R2       Adj. R2 
## ------------------------------------------------------------------------
## log_ADD_FEE     1    332.891    343.172    218.679    0.72597    0.69552 
## R_B_COST        1    333.513    343.794    219.132    0.72178    0.69087 
## OUT_STAT        1    334.111    344.392    219.570    0.71770    0.68633 
## SF_RATIO        1    341.927    352.208    225.394    0.65840    0.62045 
## GRADUAT         1    358.480    368.761    238.448    0.48849    0.43166 
## ------------------------------------------------------------------------
## 
## Step     => 3 
## Removed  => log_ADD_FEE 
## Model    => NEW10 ~ R_B_COST + SF_RATIO + GRADUAT + OUT_STAT 
## AIC      => 332.891 
## 
##                   Table: Removing Existing Variables                   
## ----------------------------------------------------------------------
## Predictor    DF      AIC        SBC       SBIC        R2       Adj. R2 
## ----------------------------------------------------------------------
## R_B_COST      1    331.621    340.189    216.806    0.72105    0.69843 
## OUT_STAT      1    332.295    340.863    217.339    0.71642    0.69343 
## SF_RATIO      1    339.927    348.495    223.430    0.65840    0.63071 
## GRADUAT       1    358.030    366.598    238.458    0.46878    0.42571 
## ----------------------------------------------------------------------
## 
## Step     => 4 
## Removed  => R_B_COST 
## Model    => NEW10 ~ SF_RATIO + GRADUAT + OUT_STAT 
## AIC      => 331.6207 
## 
##                   Table: Removing Existing Variables                   
## ----------------------------------------------------------------------
## Predictor    DF      AIC        SBC       SBIC        R2       Adj. R2 
## ----------------------------------------------------------------------
## OUT_STAT      1    330.381    337.236    215.046    0.71582    0.70087 
## SF_RATIO      1    338.121    344.975    221.632    0.65678    0.63872 
## GRADUAT       1    358.218    365.073    239.135    0.43966    0.41017 
## ----------------------------------------------------------------------
## 
## Step     => 5 
## Removed  => OUT_STAT 
## Model    => NEW10 ~ SF_RATIO + GRADUAT 
## AIC      => 330.3813 
## 
##                   Table: Removing Existing Variables                   
## ----------------------------------------------------------------------
## Predictor    DF      AIC        SBC       SBIC        R2       Adj. R2 
## ----------------------------------------------------------------------
## SF_RATIO      1    339.980    345.121    223.161    0.62291    0.61324 
## GRADUAT       1    365.570    370.711    246.575    0.29611    0.27806 
## ----------------------------------------------------------------------
## 
## 
## No more variables to be removed.
## 
## Variables Removed: 
## 
## => BOOK 
## => PH_D 
## => log_ADD_FEE 
## => R_B_COST 
## => OUT_STAT
backward_aic
## 
## 
##                               Stepwise Summary                              
## --------------------------------------------------------------------------
## Step    Variable         AIC        SBC       SBIC        R2       Adj. R2 
## --------------------------------------------------------------------------
##  0      Full Model     338.284    353.706    225.692    0.73000    0.67272 
##  1      BOOK           336.318    350.026    222.763    0.72977    0.68209 
##  2      PH_D           334.535    346.530    220.181    0.72834    0.68953 
##  3      log_ADD_FEE    332.891    343.172    217.888    0.72597    0.69552 
##  4      R_B_COST       331.621    340.189    216.109    0.72105    0.69843 
##  5      OUT_STAT       330.381    337.236    214.490    0.71582    0.70087 
## --------------------------------------------------------------------------
## 
## Final Model Output 
## ------------------
## 
##                          Model Summary                           
## ----------------------------------------------------------------
## R                        0.846       RMSE                12.336 
## R-Squared                0.716       MSE                152.183 
## Adj. R-Squared           0.701       Coef. Var           23.175 
## Pred R-Squared           0.672       AIC                330.381 
## MAE                     10.318       SBC                337.236 
## ----------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
##  AIC: Akaike Information Criteria 
##  SBC: Schwarz Bayesian Criteria 
## 
##                                 ANOVA                                 
## ---------------------------------------------------------------------
##                  Sum of                                              
##                 Squares        DF    Mean Square      F         Sig. 
## ---------------------------------------------------------------------
## Regression    15717.001         2       7858.501     47.86    0.0000 
## Residual       6239.487        38        164.197                     
## Total         21956.488        40                                    
## ---------------------------------------------------------------------
## 
##                                    Parameter Estimates                                     
## ------------------------------------------------------------------------------------------
##       model       Beta    Std. Error    Std. Beta      t        Sig       lower     upper 
## ------------------------------------------------------------------------------------------
## (Intercept)    -18.639        14.637                 -1.273    0.211    -48.270    10.992 
##    SF_RATIO     -1.760         0.499       -0.322    -3.525    0.001     -2.770    -0.749 
##     GRADUAT      1.140         0.152        0.685     7.492    0.000      0.832     1.448 
## ------------------------------------------------------------------------------------------

4.7.5.2 Forward

forward_aic<-ols_step_forward_aic(model_pr_full, progress=TRUE, details=TRUE)
## Forward Selection Method 
## ------------------------
## 
## Candidate Terms: 
## 
## 1. BOOK 
## 2. log_ADD_FEE 
## 3. R_B_COST 
## 4. PH_D 
## 5. SF_RATIO 
## 6. GRADUAT 
## 7. OUT_STAT 
## 
## 
## Step     => 0 
## Model    => NEW10 ~ 1 
## AIC      => 377.966 
## 
## Initiating stepwise selection... 
## 
##                        Table: Adding New Variables                        
## -------------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC        R2       Adj. R2  
## -------------------------------------------------------------------------
## GRADUAT         1    339.980    345.121    223.161    0.62291     0.61324 
## PH_D            1    358.045    363.186    239.617    0.41413     0.39910 
## OUT_STAT        1    362.074    367.215    243.334    0.35363     0.33706 
## SF_RATIO        1    365.570    370.711    246.575    0.29611     0.27806 
## R_B_COST        1    378.009    383.150    258.223    0.04661     0.02216 
## BOOK            1    379.594    384.735    259.720    0.00903    -0.01638 
## log_ADD_FEE     1    379.787    384.928    259.903    0.00435    -0.02118 
## -------------------------------------------------------------------------
## 
## Step     => 1 
## Added    => GRADUAT 
## Model    => NEW10 ~ GRADUAT 
## AIC      => 339.9804 
## 
##                       Table: Adding New Variables                        
## ------------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC        R2       Adj. R2 
## ------------------------------------------------------------------------
## SF_RATIO        1    330.381    337.236    215.046    0.71582    0.70087 
## OUT_STAT        1    338.121    344.975    221.632    0.65678    0.63872 
## PH_D            1    340.187    347.041    223.400    0.63905    0.62005 
## R_B_COST        1    341.456    348.310    224.490    0.62770    0.60810 
## log_ADD_FEE     1    341.939    348.793    224.905    0.62329    0.60346 
## BOOK            1    341.973    348.827    224.934    0.62297    0.60313 
## ------------------------------------------------------------------------
## 
## Step     => 2 
## Added    => SF_RATIO 
## Model    => NEW10 ~ GRADUAT + SF_RATIO 
## AIC      => 330.3813 
## 
##                       Table: Adding New Variables                        
## ------------------------------------------------------------------------
## Predictor      DF      AIC        SBC       SBIC        R2       Adj. R2 
## ------------------------------------------------------------------------
## OUT_STAT        1    331.621    340.189    216.806    0.72105    0.69843 
## log_ADD_FEE     1    332.265    340.833    217.315    0.71663    0.69366 
## PH_D            1    332.272    340.840    217.320    0.71658    0.69360 
## R_B_COST        1    332.295    340.863    217.339    0.71642    0.69343 
## BOOK            1    332.321    340.889    217.359    0.71624    0.69323 
## ------------------------------------------------------------------------
## 
## 
## No more variables to be added.
## 
## Variables Selected: 
## 
## => GRADUAT 
## => SF_RATIO
forward_aic
## 
## 
##                              Stepwise Summary                              
## -------------------------------------------------------------------------
## Step    Variable        AIC        SBC       SBIC        R2       Adj. R2 
## -------------------------------------------------------------------------
##  0      Base Model    377.966    381.393    259.401    0.00000    0.00000 
##  1      GRADUAT       339.980    345.121    223.161    0.62291    0.61324 
##  2      SF_RATIO      330.381    337.236    215.046    0.71582    0.70087 
## -------------------------------------------------------------------------
## 
## Final Model Output 
## ------------------
## 
##                          Model Summary                           
## ----------------------------------------------------------------
## R                        0.846       RMSE                12.336 
## R-Squared                0.716       MSE                152.183 
## Adj. R-Squared           0.701       Coef. Var           23.175 
## Pred R-Squared           0.672       AIC                330.381 
## MAE                     10.318       SBC                337.236 
## ----------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
##  AIC: Akaike Information Criteria 
##  SBC: Schwarz Bayesian Criteria 
## 
##                                 ANOVA                                 
## ---------------------------------------------------------------------
##                  Sum of                                              
##                 Squares        DF    Mean Square      F         Sig. 
## ---------------------------------------------------------------------
## Regression    15717.001         2       7858.501     47.86    0.0000 
## Residual       6239.487        38        164.197                     
## Total         21956.488        40                                    
## ---------------------------------------------------------------------
## 
##                                    Parameter Estimates                                     
## ------------------------------------------------------------------------------------------
##       model       Beta    Std. Error    Std. Beta      t        Sig       lower     upper 
## ------------------------------------------------------------------------------------------
## (Intercept)    -18.639        14.637                 -1.273    0.211    -48.270    10.992 
##     GRADUAT      1.140         0.152        0.685     7.492    0.000      0.832     1.448 
##    SF_RATIO     -1.760         0.499       -0.322    -3.525    0.001     -2.770    -0.749 
## ------------------------------------------------------------------------------------------

Комментарий: строим Backward и Forward AIC. Backward AIC исключил всё, кроме SF_RATIO и GRADUAT. Forward AIC включил только GRADUAT и SF_RATIO. Оба пошаговых алгоритма дали один и тот же результат.

Интерпретация параметров: число лучших 10% студентов NEW10 увеличивается с увеличением выпускников GRADUAT и уменьшается с числом студентов на преподавателя SF_RATIO.

Здесь можно сравнить по R^2, например.

4.7.6 Построение и анализ графиков

model_pr_full<-lm(NEW10 ~ GRADUAT + SF_RATIO, data=pr_data)
summary(model_pr_full)
## 
## Call:
## lm(formula = NEW10 ~ GRADUAT + SF_RATIO, data = pr_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.763  -8.954  -2.873   9.404  52.149 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -5.1676    16.4411  -0.314  0.75465    
## GRADUAT       0.9740     0.1626   5.990  2.6e-07 ***
## SF_RATIO     -1.5182     0.5621  -2.701  0.00952 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.35 on 48 degrees of freedom
##   (2 пропущенных наблюдений удалены)
## Multiple R-squared:  0.6133, Adjusted R-squared:  0.5972 
## F-statistic: 38.07 on 2 and 48 DF,  p-value: 1.247e-10
plot(model_pr_full, 1:5)

У 51-го сильное плечо (ещё видно, как линия идёт к ней)

4.7.7 Поиск outliers по Residuals vs. Deleted Residuals

plot(residuals(model_pr_full), rstudent(model_pr_full))

Комментарий: замечаем точку (51), которая сильно влияет на линейную модель.

influence.measures(model_pr_full)
## Influence measures of
##   lm(formula = NEW10 ~ GRADUAT + SF_RATIO, data = pr_data) :
## 
##       dfb.1_  dfb.GRAD  dfb.SF_R   dffit cov.r   cook.d    hat inf
## 1   0.156252 -0.045588 -2.35e-01  0.4007 0.816 4.95e-02 0.0317    
## 2  -0.142405  0.171118  7.76e-02  0.2199 1.057 1.61e-02 0.0498    
## 3   0.007595 -0.007708 -9.15e-05  0.0145 1.099 7.16e-05 0.0308    
## 4   0.007442  0.011854 -7.56e-02 -0.1148 1.115 4.46e-03 0.0579    
## 5  -0.062742  0.118154 -4.21e-02  0.2030 1.079 1.38e-02 0.0554    
## 6   0.209642 -0.134914 -3.80e-01 -0.4260 1.077 5.98e-02 0.1036    
## 7  -0.010436  0.018846 -3.62e-02 -0.0813 1.093 2.24e-03 0.0360    
## 8  -0.156587  0.114409  8.70e-02 -0.3109 0.849 3.03e-02 0.0231    
## 9  -0.157795  0.122821  1.25e-01 -0.2056 1.037 1.41e-02 0.0387    
## 10  0.025289 -0.043555  6.54e-03 -0.0741 1.102 1.86e-03 0.0411    
## 12 -0.031867  0.027892  2.43e-02 -0.0335 1.206 3.81e-04 0.1170   *
## 13  0.007036 -0.029448  3.02e-02 -0.0739 1.117 1.85e-03 0.0523    
## 14  0.148220 -0.160722 -2.28e-02  0.2043 1.101 1.40e-02 0.0676    
## 15 -0.440485  0.369605  3.53e-01 -0.4841 0.980 7.54e-02 0.0802    
## 16 -0.020974  0.058786 -3.49e-02  0.1383 1.069 6.44e-03 0.0357    
## 17 -0.013268 -0.014281  4.90e-02 -0.0892 1.099 2.70e-03 0.0421    
## 18 -0.165762  0.166156  1.42e-01  0.2044 1.139 1.41e-02 0.0899    
## 19 -0.037008  0.028235  2.21e-02 -0.0639 1.082 1.39e-03 0.0251    
## 20 -0.014143  0.000378  3.18e-02 -0.0433 1.136 6.37e-04 0.0640    
## 21 -0.169351  0.147972  1.93e-01  0.2218 1.183 1.66e-02 0.1200    
## 22  0.024521 -0.044146 -1.97e-02 -0.1356 1.034 6.15e-03 0.0219    
## 24 -0.233928  0.276537  1.44e-01  0.3683 0.934 4.35e-02 0.0451    
## 25 -0.081062  0.088996  1.72e-02 -0.1024 1.203 3.56e-03 0.1191   *
## 26  0.071188 -0.084065 -4.77e-02 -0.1185 1.086 4.74e-03 0.0399    
## 27 -0.182509  0.145727  1.67e-01 -0.2044 1.152 1.41e-02 0.0979    
## 28 -0.067432 -0.008354  1.65e-01 -0.2386 1.063 1.90e-02 0.0563    
## 29 -0.008973  0.028263 -2.49e-02  0.0639 1.127 1.39e-03 0.0585    
## 30 -0.110949  0.144068  3.66e-02  0.1853 1.086 1.15e-02 0.0551    
## 31  0.084879 -0.059446 -1.54e-01 -0.1872 1.099 1.18e-02 0.0624    
## 32  0.022119 -0.040726  1.33e-02 -0.0689 1.121 1.61e-03 0.0550    
## 33 -0.016792  0.037689 -1.59e-02  0.0784 1.095 2.09e-03 0.0372    
## 34  0.154239 -0.125010 -2.45e-01 -0.3179 0.989 3.30e-02 0.0484    
## 35 -0.052304  0.065696 -2.44e-02 -0.1143 1.110 4.43e-03 0.0545    
## 36  0.189408 -0.151498 -1.43e-01  0.2417 1.014 1.93e-02 0.0395    
## 37 -0.029891 -0.001702  1.42e-01  0.2076 1.068 1.44e-02 0.0514    
## 38 -0.105057  0.106146  9.90e-03 -0.1824 1.035 1.11e-02 0.0329    
## 39 -0.060391  0.032813  7.43e-02 -0.1045 1.091 3.70e-03 0.0397    
## 40 -0.018011  0.057624 -5.00e-02  0.1320 1.105 5.89e-03 0.0543    
## 41  0.373278 -0.291574 -3.66e-01  0.4214 1.143 5.90e-02 0.1335    
## 42  0.041672 -0.031904 -2.39e-02  0.0735 1.077 1.83e-03 0.0247    
## 43 -0.091131  0.106561 -1.98e-02 -0.1722 1.081 9.97e-03 0.0494    
## 44  0.288062 -0.292022 -2.68e-01 -0.4058 0.948 5.29e-02 0.0556    
## 45 -0.015808  0.050662 -3.76e-02  0.1224 1.083 5.06e-03 0.0394    
## 46 -0.062436  0.092305  4.33e-03  0.1352 1.091 6.17e-03 0.0466    
## 47  0.024276 -0.009007 -3.60e-02  0.0538 1.102 9.82e-04 0.0377    
## 48 -0.073400  0.113328  1.46e-02  0.1985 1.019 1.31e-02 0.0315    
## 49  0.000713  0.003235 -2.22e-02 -0.0423 1.099 6.09e-04 0.0335    
## 50 -0.004980  0.008760 -1.19e-02 -0.0260 1.121 2.29e-04 0.0506    
## 51  1.078464 -1.600184  8.12e-01  2.5174 0.450 1.48e+00 0.2292   *
## 52  0.047770 -0.034990 -1.09e-01 -0.1676 1.051 9.41e-03 0.0348    
## 53  0.027573 -0.015210 -5.76e-02 -0.0660 1.184 1.48e-03 0.1022
row.names(model_pr_full$model)
##  [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "12" "13" "14" "15" "16"
## [16] "17" "18" "19" "20" "21" "22" "24" "25" "26" "27" "28" "29" "30" "31" "32"
## [31] "33" "34" "35" "36" "37" "38" "39" "40" "41" "42" "43" "44" "45" "46" "47"
## [46] "48" "49" "50" "51" "52" "53"
outlier_cook<-pr_data|>dplyr::filter(NEW10==48&SF_RATIO==20.5)
outlier_cook

Комментарий: выбросом оказался Brigham Young University.

4.7.8 Поиск outliers по Cook и Mahalanobis

4.7.8.1 Cook

# Рассчитаем Cook's distance для модели с логарифмами
cooks_vals <- cooks.distance(model_pr_full)
cook_df <- data.frame(
  obs  = row.names(model_pr_full$model),
  cook = cooks_vals
) %>%
  arrange(desc(cook))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(cook_df))


plot(
  x, cook_df$cook,
  type  = "h",
  main  = "Cook's Distance",
  xlab  = "Номер наблюдения (из исходного df)",
  ylab  = "Cook's distance",
  xaxt  = "n"                # убираем автоматическую ось X
)
axis(side = 1, at = x, labels = cook_df$obs, las = 2, cex.axis = 0.8)
abline(h = 1, col = "red", lty = 2)

Комментарий: высокие значения у 29 и 38.

4.7.8.2 Mahalanobis

# 7.2. Расстояние Махаланобиса (Mahalanobis Distance)
# Расстояние Махаланобиса в пространстве предикторов (регрессоров)
mahalanobis_distance <- mahalanobis(model.matrix(model_pr_full)[,-1],
                                     center = colMeans(model.matrix(model_pr_full)[,-1]),
                                     cov = cov(model.matrix(model_pr_full)[,-1]))
mah_df <- data.frame(
  obs  = row.names(model_pr_full$model),
  mah = mahalanobis_distance
) %>%
  arrange(desc(mah))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(mah_df))
alpha<-0.05

# Scatter plot для расстояния Махаланобиса
plot(
  x, mah_df$mah,
  type  = "h",
  main  = "Mahalanobis's Distance",
  xlab  = "Номер наблюдения (из исходного df)",
  ylab  = "Mahalanobis's distance",
  xaxt  = "n"                
)
axis(side = 1, at = x, labels = mah_df$obs, las = 2, cex.axis = 0.8)
abline(h = qchisq(p = 1 - alpha, df = nrow(mah_df)), col = "red", lty = 2)

Комментарий: 29-ое наблюдение имеет большое расстояние по Махаланобису.

4.7.9 Линейная модель без выброса

pr_data_outlier<-pr_data[-c(51),]
plot_numeric_pairs(pr_data_outlier)

model_pr_outlier<-lm(NEW10~GRADUAT + SF_RATIO, data=pr_data_outlier)
summary(model_pr_full|>lm.beta())
## 
## Call:
## lm(formula = NEW10 ~ GRADUAT + SF_RATIO, data = pr_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.763  -8.954  -2.873   9.404  52.149 
## 
## Coefficients:
##             Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept)  -5.1676           NA    16.4411  -0.314  0.75465    
## GRADUAT       0.9740       0.6124     0.1626   5.990  2.6e-07 ***
## SF_RATIO     -1.5182      -0.2762     0.5621  -2.701  0.00952 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.35 on 48 degrees of freedom
##   (2 пропущенных наблюдений удалены)
## Multiple R-squared:  0.6133, Adjusted R-squared:  0.5972 
## F-statistic: 38.07 on 2 and 48 DF,  p-value: 1.247e-10
summary(model_pr_outlier|>lm.beta())
## 
## Call:
## lm(formula = NEW10 ~ GRADUAT + SF_RATIO, data = pr_data_outlier)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -27.6995  -8.8741   0.3934   8.0222  31.4411 
## 
## Coefficients:
##             Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept) -20.0312           NA    14.1534  -1.415 0.163576    
## GRADUAT       1.1921       0.6758     0.1443   8.263 1.03e-10 ***
## SF_RATIO     -1.9009      -0.3250     0.4784  -3.973 0.000242 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.87 on 47 degrees of freedom
##   (2 пропущенных наблюдений удалены)
## Multiple R-squared:  0.7332, Adjusted R-squared:  0.7218 
## F-statistic: 64.57 on 2 and 47 DF,  p-value: 3.284e-14

Комментарий: замечаем существенную разницу в R^2, p-value и значимости SF_RATIO.

plot(model_pr_outlier, 1:5)

AIC(model_pr_full)
## [1] 428.221
AIC(model_pr_outlier)
## [1] 402.2787

Можем заметить, что модель без выброса лучше по AIC, чем с выбросом.

4.7.10 Итог: прогнозирование

set.seed(456) # Другое зерно для генерации новых значений
max_GRADUAT<-max(pr_data$GRADUAT, na.rm=TRUE)
max_SF_RATIO<-max(pr_data$SF_RATIO, na.rm=TRUE)
new_data <- data.frame(
  GRADUAT = max_GRADUAT - 1:8,
  SF_RATIO = max_SF_RATIO - 1:8
)

# Прогнозируем значение NEW10
predicted_values <- predict(model_pr_outlier, newdata = new_data)
cat("Прогнозируемые значения NEW10:\n")
## Прогнозируемые значения NEW10:
print(predicted_values)
##        1        2        3        4        5        6        7        8 
## 59.72933 60.43809 61.14684 61.85560 62.56435 63.27311 63.98186 64.69061
# 6. Доверительные и предсказательные интервалы

# Доверительный интервал (confidence interval) для среднего значения NEW10
confidence_intervals <- predict(model_pr_outlier, newdata = new_data, interval = "confidence", level = 0.95)
print("Доверительные интервалы (95%):\n")
## [1] "Доверительные интервалы (95%):\n"
print(confidence_intervals)
##        fit      lwr      upr
## 1 59.72933 46.77615 72.68252
## 2 60.43809 48.53109 72.34508
## 3 61.14684 50.27377 72.01991
## 4 61.85560 52.00033 71.71086
## 5 62.56435 53.70520 71.42350
## 6 63.27311 55.38017 71.16604
## 7 63.98186 57.01279 70.95093
## 8 64.69061 58.58381 70.79742
# Предсказательный интервал (prediction interval) для отдельного наблюдения NEW10
prediction_intervals <- predict(model_pr_outlier, newdata = new_data, interval = "prediction", level = 0.95)
print("Предсказательные интервалы (95%):\n")
## [1] "Предсказательные интервалы (95%):\n"
print(prediction_intervals)
##        fit      lwr      upr
## 1 59.72933 30.78140 88.67727
## 2 60.43809 31.94292 88.93325
## 3 61.14684 33.06800 89.22569
## 4 61.85560 34.15498 89.55621
## 5 62.56435 35.20229 89.92641
## 6 63.27311 36.20843 90.33778
## 7 63.98186 37.17205 90.79167
## 8 64.69061 38.09191 91.28932
# Создаем data.frame для прогнозов и интервалов
predictions_df <- data.frame(
  GRADUAT = new_data$GRADUAT,
  SF_RATIO = new_data$SF_RATIO,
  Predicted = predicted_values,
  Lower_Conf = confidence_intervals[, "lwr"],
  Upper_Conf = confidence_intervals[, "upr"],
  Lower_Pred = prediction_intervals[, "lwr"],
  Upper_Pred = prediction_intervals[, "upr"]
)
predictions_df
# Визуализация (только для демонстрации, нужен более сложный график для нескольких предикторов)
ggplot(predictions_df, aes(x = 1:nrow(predictions_df), y = Predicted)) + # Условная ось x, т.к. несколько предикторов
  geom_point() +
  geom_errorbar(aes(ymin = Lower_Conf, ymax = Upper_Conf), color = "blue", width = 0.2) + # Доверительные интервалы
  geom_errorbar(aes(ymin = Lower_Pred, ymax = Upper_Pred), color = "red", width = 0.2) + # Предсказательные интервалы
  labs(title = "Прогнозы NEW10 с доверительными и предсказательными интервалами",
       x = "Прогноз (номер)",
       y = "NEW10") +
  scale_x_continuous(breaks = 1:nrow(predictions_df)) + # Подписи для оси x
  theme_bw()

  # annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.9, label = "Синие: Доверительные интервалы", color = "blue", hjust = 0) +
  # annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.8, label = "Красные: Предсказательные интервалы", color = "red", hjust = 0)

4.8 Анализ общественных вузов

pub_data<-col_I_regr_split$public
pub_data<-pub_data|>mutate(NEW10 = log(NEW10))|>rename(log_NEW10 = NEW10)
print(paste0("Число наблюдений: ", dim(pub_data)[1]))
## [1] "Число наблюдений: 123"
colSums(is.na(pub_data))
##        ...1 log_ADD_FEE        BOOK    R_B_COST    OUT_STAT   log_NEW10 
##           0          33           2           0           1          15 
##    SF_RATIO        PH_D     GRADUAT       PPIND 
##           0           5           4           0

Много пропусков по log_ADD_FEE (около четверти).

4.8.1 График

plot_numeric_pairs(pub_data, "Общественные вузы")

4.8.2 Выбросы “на глаз”

outlier_pub_out_eye<-pub_data[pub_data$OUT_STAT>14000,]
outlier_pub_out_eye
plot_numeric_pairs(pub_data, "Общественные вузы", list(red = outlier_pub_out_eye))

4.8.3 Линейная модель

library(lm.beta)
model_pub_full <- lm(log_NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT, 
                    data = pub_data)
summary(model_pub_full|>lm.beta())
## 
## Call:
## lm(formula = log_NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + 
##     SF_RATIO + GRADUAT + OUT_STAT, data = pub_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.23717 -0.22303  0.04705  0.22529  0.89431 
## 
## Coefficients:
##               Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept) -4.821e-02           NA  8.271e-01  -0.058   0.9537    
## BOOK         8.739e-04    1.632e-01  5.489e-04   1.592   0.1161    
## log_ADD_FEE  4.078e-03    7.717e-03  5.368e-02   0.076   0.9397    
## R_B_COST     8.749e-05    1.289e-01  7.275e-05   1.203   0.2334    
## PH_D         2.303e-02    2.706e-01  8.974e-03   2.566   0.0126 *  
## SF_RATIO    -9.484e-03   -7.112e-02  1.289e-02  -0.736   0.4645    
## GRADUAT      1.970e-02    5.464e-01  4.132e-03   4.768 1.07e-05 ***
## OUT_STAT    -5.986e-05   -2.567e-01  2.867e-05  -2.088   0.0407 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4196 on 66 degrees of freedom
##   (49 пропущенных наблюдений удалены)
## Multiple R-squared:  0.4814, Adjusted R-squared:  0.4264 
## F-statistic: 8.752 on 7 and 66 DF,  p-value: 1.449e-07

4.8.4 Доверительные эллипсы

plot_confidence_ellipse(model_pub_full, c(2, 3))

plot_confidence_ellipse(model_pub_full, c(4, 5))

plot_confidence_ellipse(model_pub_full, c(6, 7))

plot_confidence_ellipse(model_pub_full, c(8, 2))

# Автор: -Евгений-

# Вычисление множественной корреляции
calc_multicollinearity <- function(var, data) {
  X <- as.matrix(data[, setdiff(names(data), var)])  # Матрица остальных регрессоров
  y <- as.matrix(data[[var]])  # Текущий регрессор как зависимая переменная
  model <- lm(y ~ X) 
  sqrt(summary(model)$r.squared)  
}

calc_partial_correlation <- function(var, data) {
  # Остатки регрессии текущего регрессора на остальные регрессоры
  data<-data|>na.omit()
  model_X <- lm(data[[var]] ~ ., data = data[, !names(data) %in% c(var, "log_NEW10")])
  res_X <- residuals(model_X)
  
  # Остатки регрессии Y на остальные регрессоры
  model_Y <- lm(log_NEW10 ~ ., data = data[, !names(data) %in% var])
  res_Y <- residuals(model_Y)
  
  cor(res_X, res_Y, use = "complete.obs")  # Корреляция остатков
}

regressors <- setdiff(names(pub_data), c("...1", "PPIND", "log_NEW10"))
regressors_full<-c(regressors, "log_NEW10")
redundancy_table <- data.frame(
  Регрессор = regressors,
  Множественная_корреляция = sapply(regressors, calc_multicollinearity, pub_data[, regressors_full]),
  Частная_корреляция = sapply(regressors, calc_partial_correlation, pub_data[,regressors_full])
)

print(redundancy_table)
##               Регрессор Множественная_корреляция Частная_корреляция
## log_ADD_FEE log_ADD_FEE                0.4883582        0.009351778
## BOOK               BOOK                0.5291391        0.192313348
## R_B_COST       R_B_COST                0.5747925        0.146434350
## OUT_STAT       OUT_STAT                0.7158033       -0.248923976
## SF_RATIO       SF_RATIO                0.4073614       -0.090198497
## PH_D               PH_D                0.5981262        0.301176912
## GRADUAT         GRADUAT                0.7449868        0.506130402

Исключим log_ADD_FEE и SF_RATIO (так как у них низкие частные корреляции)

AIC(model_pub_full)
## [1] 90.99499
summary(model_pub_full)
## 
## Call:
## lm(formula = log_NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + 
##     SF_RATIO + GRADUAT + OUT_STAT, data = pub_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.23717 -0.22303  0.04705  0.22529  0.89431 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.821e-02  8.271e-01  -0.058   0.9537    
## BOOK         8.739e-04  5.489e-04   1.592   0.1161    
## log_ADD_FEE  4.078e-03  5.368e-02   0.076   0.9397    
## R_B_COST     8.749e-05  7.275e-05   1.203   0.2334    
## PH_D         2.303e-02  8.974e-03   2.566   0.0126 *  
## SF_RATIO    -9.484e-03  1.289e-02  -0.736   0.4645    
## GRADUAT      1.970e-02  4.132e-03   4.768 1.07e-05 ***
## OUT_STAT    -5.986e-05  2.867e-05  -2.088   0.0407 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4196 on 66 degrees of freedom
##   (49 пропущенных наблюдений удалены)
## Multiple R-squared:  0.4814, Adjusted R-squared:  0.4264 
## F-statistic: 8.752 on 7 and 66 DF,  p-value: 1.449e-07
model_pub_nono<-lm(log_NEW10 ~ R_B_COST + PH_D + GRADUAT + OUT_STAT + BOOK, data=pub_data)
summary(model_pub_full|>lm.beta())
## 
## Call:
## lm(formula = log_NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + 
##     SF_RATIO + GRADUAT + OUT_STAT, data = pub_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.23717 -0.22303  0.04705  0.22529  0.89431 
## 
## Coefficients:
##               Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept) -4.821e-02           NA  8.271e-01  -0.058   0.9537    
## BOOK         8.739e-04    1.632e-01  5.489e-04   1.592   0.1161    
## log_ADD_FEE  4.078e-03    7.717e-03  5.368e-02   0.076   0.9397    
## R_B_COST     8.749e-05    1.289e-01  7.275e-05   1.203   0.2334    
## PH_D         2.303e-02    2.706e-01  8.974e-03   2.566   0.0126 *  
## SF_RATIO    -9.484e-03   -7.112e-02  1.289e-02  -0.736   0.4645    
## GRADUAT      1.970e-02    5.464e-01  4.132e-03   4.768 1.07e-05 ***
## OUT_STAT    -5.986e-05   -2.567e-01  2.867e-05  -2.088   0.0407 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4196 on 66 degrees of freedom
##   (49 пропущенных наблюдений удалены)
## Multiple R-squared:  0.4814, Adjusted R-squared:  0.4264 
## F-statistic: 8.752 on 7 and 66 DF,  p-value: 1.449e-07
AIC(model_pub_full)
## [1] 90.99499
summary(model_pub_nono|>lm.beta())
## 
## Call:
## lm(formula = log_NEW10 ~ R_B_COST + PH_D + GRADUAT + OUT_STAT + 
##     BOOK, data = pub_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.21996 -0.24955 -0.00251  0.24842  0.81900 
## 
## Coefficients:
##               Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept)  7.390e-02           NA  5.501e-01   0.134   0.8934    
## R_B_COST     1.272e-04    1.806e-01  6.475e-05   1.964   0.0525 .  
## PH_D         1.748e-02    2.156e-01  7.462e-03   2.342   0.0213 *  
## GRADUAT      1.808e-02    4.905e-01  3.432e-03   5.268 9.03e-07 ***
## OUT_STAT    -5.027e-05   -2.003e-01  2.389e-05  -2.104   0.0381 *  
## BOOK         1.044e-03    1.903e-01  4.510e-04   2.315   0.0228 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4225 on 92 degrees of freedom
##   (25 пропущенных наблюдений удалены)
## Multiple R-squared:  0.4772, Adjusted R-squared:  0.4488 
## F-statistic:  16.8 on 5 and 92 DF,  p-value: 9.211e-12
AIC(model_pub_nono)
## [1] 117.0556

Комментарий: AIC увеличился, p-value уменьшился. Теперь каждый признак, кроме R_B_COST является значимым (при alpha=0.05).

regressors <- setdiff(names(pub_data), c("...1", "PPIND", "log_NEW10", "log_ADD_FEE", "SF_RATIO"))
regressors_full<-c(regressors, "log_NEW10")
redundancy_table <- data.frame(
  Регрессор = regressors,
  Множественная_корреляция = sapply(regressors, calc_multicollinearity, pub_data[, regressors_full]),
  Частная_корреляция = sapply(regressors, calc_partial_correlation, pub_data[,regressors_full])
)

print(redundancy_table)
##          Регрессор Множественная_корреляция Частная_корреляция
## BOOK          BOOK                0.4535409          0.2346185
## R_B_COST  R_B_COST                0.5958627          0.2006126
## OUT_STAT  OUT_STAT                0.6338016         -0.2142731
## PH_D          PH_D                0.6059169          0.2372086
## GRADUAT    GRADUAT                0.7046461          0.4813910

Тут уже не так понятно. Перейдём к пошаговому AIC.

AIC(model_pub_full)
## [1] 90.99499
AIC(model_pub_nono)
## [1] 117.0556

4.8.5 OLS. Пошаговый AIC.

4.8.5.1 Backward

library(olsrr)
backward_aic<-ols_step_backward_aic(model_pub_full, progress=TRUE, details=TRUE)
## Backward Elimination Method 
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. BOOK 
## 2. log_ADD_FEE 
## 3. R_B_COST 
## 4. PH_D 
## 5. SF_RATIO 
## 6. GRADUAT 
## 7. OUT_STAT 
## 
## 
## Step     => 0 
## Model    => log_NEW10 ~ BOOK + log_ADD_FEE + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT 
## AIC      => 90.99499 
## 
## Initiating stepwise selection... 
## 
##                    Table: Removing Existing Variables                     
## -------------------------------------------------------------------------
## Predictor      DF      AIC        SBC        SBIC        R2       Adj. R2 
## -------------------------------------------------------------------------
## log_ADD_FEE     1     89.001    107.434    -119.335    0.48135    0.43490 
## SF_RATIO        1     89.599    108.032    -118.859    0.47714    0.43032 
## R_B_COST        1     90.599    109.032    -118.062    0.47003    0.42257 
## BOOK            1     91.784    110.216    -117.115    0.46148    0.41325 
## OUT_STAT        1     93.728    112.161    -115.555    0.44714    0.39763 
## PH_D            1     96.032    114.464    -113.699    0.42966    0.37858 
## GRADUAT         1    110.895    129.327    -101.488    0.30279    0.24035 
## -------------------------------------------------------------------------
## 
## Step     => 1 
## Removed  => log_ADD_FEE 
## Model    => log_NEW10 ~ BOOK + R_B_COST + PH_D + SF_RATIO + GRADUAT + OUT_STAT 
## AIC      => 89.00146 
## 
##                   Table: Removing Existing Variables                    
## -----------------------------------------------------------------------
## Predictor    DF      AIC        SBC        SBIC        R2       Adj. R2 
## -----------------------------------------------------------------------
## SF_RATIO      1     87.604    103.733    -121.080    0.47711    0.43866 
## R_B_COST      1     88.775    104.904    -120.111    0.46877    0.42970 
## BOOK          1     89.835    105.964    -119.232    0.46110    0.42148 
## OUT_STAT      1     91.997    108.125    -117.436    0.44513    0.40433 
## PH_D          1     94.185    110.314    -115.609    0.42847    0.38645 
## GRADUAT       1    111.437    127.565    -100.972    0.27842    0.22537 
## -----------------------------------------------------------------------
## 
## Step     => 2 
## Removed  => SF_RATIO 
## Model    => log_NEW10 ~ BOOK + R_B_COST + PH_D + GRADUAT + OUT_STAT 
## AIC      => 87.60406 
## 
##                   Table: Removing Existing Variables                    
## -----------------------------------------------------------------------
## Predictor    DF      AIC        SBC        SBIC        R2       Adj. R2 
## -----------------------------------------------------------------------
## R_B_COST      1     87.268    101.092    -121.878    0.46522    0.43422 
## BOOK          1     88.652    102.476    -120.689    0.45512    0.42354 
## OUT_STAT      1     90.007    103.831    -119.523    0.44505    0.41288 
## PH_D          1     94.291    108.115    -115.824    0.41198    0.37789 
## GRADUAT       1    109.486    123.311    -102.540    0.27794    0.23608 
## -----------------------------------------------------------------------
## 
## Step     => 3 
## Removed  => R_B_COST 
## Model    => log_NEW10 ~ BOOK + PH_D + GRADUAT + OUT_STAT 
## AIC      => 87.26753 
## 
##                   Table: Removing Existing Variables                    
## -----------------------------------------------------------------------
## Predictor    DF      AIC        SBC        SBIC        R2       Adj. R2 
## -----------------------------------------------------------------------
## OUT_STAT      1     88.405     99.925    -121.264    0.44206    0.41815 
## BOOK          1     89.065    100.585    -120.677    0.43706    0.41294 
## PH_D          1     96.921    108.442    -113.661    0.37401    0.34718 
## GRADUAT       1    109.030    120.550    -102.753    0.26272    0.23112 
## -----------------------------------------------------------------------
## 
## 
## No more variables to be removed.
## 
## Variables Removed: 
## 
## => log_ADD_FEE 
## => SF_RATIO 
## => R_B_COST
backward_aic
## 
## 
##                               Stepwise Summary                              
## --------------------------------------------------------------------------
## Step    Variable        AIC        SBC        SBIC        R2       Adj. R2 
## --------------------------------------------------------------------------
##  0      Full Model     90.995    111.732    -117.098    0.48139    0.42639 
##  1      log_ADD_FEE    89.001    107.434    -119.561    0.48135    0.43490 
##  2      SF_RATIO       87.604    103.733    -121.356    0.47711    0.43866 
##  3      R_B_COST       87.268    101.092    -122.021    0.46522    0.43422 
## --------------------------------------------------------------------------
## 
## Final Model Output 
## ------------------
## 
##                          Model Summary                          
## ---------------------------------------------------------------
## R                       0.682       RMSE                 0.402 
## R-Squared               0.465       MSE                  0.162 
## Adj. R-Squared          0.434       Coef. Var           12.482 
## Pred R-Squared          0.385       AIC                 87.268 
## MAE                     0.300       SBC                101.092 
## ---------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
##  AIC: Akaike Information Criteria 
##  SBC: Schwarz Bayesian Criteria 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression     10.423         4          2.606    15.006    0.0000 
## Residual       11.981        69          0.174                     
## Total          22.403        73                                    
## -------------------------------------------------------------------
## 
##                                   Parameter Estimates                                   
## ---------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower    upper 
## ---------------------------------------------------------------------------------------
## (Intercept)    -0.435         0.670                 -0.649    0.519    -1.771    0.902 
##        BOOK     0.001         0.001        0.188     1.906    0.061     0.000    0.002 
##        PH_D     0.028         0.008        0.329     3.431    0.001     0.012    0.044 
##     GRADUAT     0.019         0.004        0.538     5.112    0.000     0.012    0.027 
##    OUT_STAT     0.000         0.000       -0.187    -1.729    0.088     0.000    0.000 
## ---------------------------------------------------------------------------------------

4.8.5.2 Forward

forward_aic<-ols_step_forward_aic(model_pub_full, progress=TRUE, details=TRUE)
## Forward Selection Method 
## ------------------------
## 
## Candidate Terms: 
## 
## 1. BOOK 
## 2. log_ADD_FEE 
## 3. R_B_COST 
## 4. PH_D 
## 5. SF_RATIO 
## 6. GRADUAT 
## 7. OUT_STAT 
## 
## 
## Step     => 0 
## Model    => log_NEW10 ~ 1 
## AIC      => 125.5841 
## 
## Initiating stepwise selection... 
## 
##                        Table: Adding New Variables                         
## --------------------------------------------------------------------------
## Predictor      DF      AIC        SBC        SBIC        R2       Adj. R2  
## --------------------------------------------------------------------------
## GRADUAT         1    107.184    114.097    -103.864    0.24094     0.23040 
## PH_D            1    107.752    114.664    -103.325    0.23509     0.22447 
## R_B_COST        1    119.631    126.543     -92.031    0.10191     0.08943 
## BOOK            1    121.310    128.222     -90.431    0.08130     0.06854 
## log_ADD_FEE     1    122.064    128.977     -89.711    0.07188     0.05899 
## SF_RATIO        1    126.507    133.419     -85.472    0.01445     0.00076 
## OUT_STAT        1    127.284    134.197     -84.730    0.00404    -0.00979 
## --------------------------------------------------------------------------
## 
## Step     => 1 
## Added    => GRADUAT 
## Model    => log_NEW10 ~ GRADUAT 
## AIC      => 107.1844 
## 
##                        Table: Adding New Variables                        
## -------------------------------------------------------------------------
## Predictor      DF      AIC        SBC        SBIC        R2       Adj. R2 
## -------------------------------------------------------------------------
## PH_D            1     92.767    101.983    -117.502    0.39197    0.37484 
## BOOK            1     96.918    106.134    -113.678    0.35689    0.33877 
## OUT_STAT        1    104.070    113.286    -107.072    0.29163    0.27168 
## R_B_COST        1    105.322    114.539    -105.913    0.27954    0.25924 
## SF_RATIO        1    107.925    117.141    -103.501    0.25375    0.23273 
## log_ADD_FEE     1    108.311    117.527    -103.142    0.24984    0.22871 
## -------------------------------------------------------------------------
## 
## Step     => 2 
## Added    => PH_D 
## Model    => log_NEW10 ~ GRADUAT + PH_D 
## AIC      => 92.76671 
## 
##                       Table: Adding New Variables                        
## ------------------------------------------------------------------------
## Predictor      DF     AIC        SBC        SBIC        R2       Adj. R2 
## ------------------------------------------------------------------------
## BOOK            1    88.405     99.925    -121.264    0.44206    0.41815 
## OUT_STAT        1    89.065    100.585    -120.677    0.43706    0.41294 
## R_B_COST        1    94.202    105.723    -116.094    0.39659    0.37073 
## log_ADD_FEE     1    94.654    106.174    -115.690    0.39290    0.36688 
## SF_RATIO        1    94.759    106.280    -115.596    0.39203    0.36597 
## ------------------------------------------------------------------------
## 
## Step     => 3 
## Added    => BOOK 
## Model    => log_NEW10 ~ GRADUAT + PH_D + BOOK 
## AIC      => 88.40477 
## 
##                       Table: Adding New Variables                        
## ------------------------------------------------------------------------
## Predictor      DF     AIC        SBC        SBIC        R2       Adj. R2 
## ------------------------------------------------------------------------
## OUT_STAT        1    87.268    101.092    -121.878    0.46522    0.43422 
## R_B_COST        1    90.007    103.831    -119.523    0.44505    0.41288 
## log_ADD_FEE     1    90.016    103.840    -119.515    0.44498    0.41281 
## SF_RATIO        1    90.387    104.212    -119.195    0.44219    0.40986 
## ------------------------------------------------------------------------
## 
## Step     => 4 
## Added    => OUT_STAT 
## Model    => log_NEW10 ~ GRADUAT + PH_D + BOOK + OUT_STAT 
## AIC      => 87.26753 
## 
##                       Table: Adding New Variables                        
## ------------------------------------------------------------------------
## Predictor      DF     AIC        SBC        SBIC        R2       Adj. R2 
## ------------------------------------------------------------------------
## R_B_COST        1    87.604    103.733    -121.080    0.47711    0.43866 
## SF_RATIO        1    88.775    104.904    -120.111    0.46877    0.42970 
## log_ADD_FEE     1    89.109    105.237    -119.834    0.46636    0.42713 
## ------------------------------------------------------------------------
## 
## 
## No more variables to be added.
## 
## Variables Selected: 
## 
## => GRADUAT 
## => PH_D 
## => BOOK 
## => OUT_STAT
forward_aic
## 
## 
##                               Stepwise Summary                              
## --------------------------------------------------------------------------
## Step    Variable        AIC        SBC        SBIC        R2       Adj. R2 
## --------------------------------------------------------------------------
##  0      Base Model    125.584    130.192     -85.606    0.00000    0.00000 
##  1      GRADUAT       107.184    114.097    -103.864    0.24094    0.23040 
##  2      PH_D           92.767    101.983    -117.502    0.39197    0.37484 
##  3      BOOK           88.405     99.925    -121.264    0.44206    0.41815 
##  4      OUT_STAT       87.268    101.092    -121.878    0.46522    0.43422 
## --------------------------------------------------------------------------
## 
## Final Model Output 
## ------------------
## 
##                          Model Summary                          
## ---------------------------------------------------------------
## R                       0.682       RMSE                 0.402 
## R-Squared               0.465       MSE                  0.162 
## Adj. R-Squared          0.434       Coef. Var           12.482 
## Pred R-Squared          0.385       AIC                 87.268 
## MAE                     0.300       SBC                101.092 
## ---------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
##  AIC: Akaike Information Criteria 
##  SBC: Schwarz Bayesian Criteria 
## 
##                                ANOVA                                
## -------------------------------------------------------------------
##                Sum of                                              
##               Squares        DF    Mean Square      F         Sig. 
## -------------------------------------------------------------------
## Regression     10.423         4          2.606    15.006    0.0000 
## Residual       11.981        69          0.174                     
## Total          22.403        73                                    
## -------------------------------------------------------------------
## 
##                                   Parameter Estimates                                   
## ---------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig      lower    upper 
## ---------------------------------------------------------------------------------------
## (Intercept)    -0.435         0.670                 -0.649    0.519    -1.771    0.902 
##     GRADUAT     0.019         0.004        0.538     5.112    0.000     0.012    0.027 
##        PH_D     0.028         0.008        0.329     3.431    0.001     0.012    0.044 
##        BOOK     0.001         0.001        0.188     1.906    0.061     0.000    0.002 
##    OUT_STAT     0.000         0.000       -0.187    -1.729    0.088     0.000    0.000 
## ---------------------------------------------------------------------------------------

Комментарий: строим Backward и Forward AIC. Backward AIC исключил всё, кроме SF_RATIO и GRADUAT. Forward AIC включил только GRADUAT и SF_RATIO. Оба пошаговых алгоритма дали один и тот же результат.

Интерпретация параметров: число лучших 10% студентов NEW10 увеличивается с увеличением выпускников GRADUAT, числом докторов (PH_D), ценой за бронирование (BOOK) и уменьшается со стоимостью оплаты OUT_STAT.

Здесь можно сравнить по R^2, например.

4.8.6 Построение и анализ графиков

model_pub_full<-lm(log_NEW10 ~ GRADUAT + PH_D + BOOK + OUT_STAT, data=pub_data)
summary(model_pub_full)
## 
## Call:
## lm(formula = log_NEW10 ~ GRADUAT + PH_D + BOOK + OUT_STAT, data = pub_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.21130 -0.24350 -0.01198  0.24097  0.92728 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.967e-02  5.542e-01  -0.108  0.91449    
## GRADUAT      1.825e-02  3.483e-03   5.241 9.95e-07 ***
## PH_D         2.271e-02  7.077e-03   3.209  0.00183 ** 
## BOOK         1.152e-03  4.545e-04   2.535  0.01290 *  
## OUT_STAT    -3.446e-05  2.284e-05  -1.509  0.13474    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4289 on 93 degrees of freedom
##   (25 пропущенных наблюдений удалены)
## Multiple R-squared:  0.4553, Adjusted R-squared:  0.4319 
## F-statistic: 19.43 on 4 and 93 DF,  p-value: 1.196e-11
AIC(model_pub_full)
## [1] 119.0812
summary(model_pub_nono)
## 
## Call:
## lm(formula = log_NEW10 ~ R_B_COST + PH_D + GRADUAT + OUT_STAT + 
##     BOOK, data = pub_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.21996 -0.24955 -0.00251  0.24842  0.81900 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.390e-02  5.501e-01   0.134   0.8934    
## R_B_COST     1.272e-04  6.475e-05   1.964   0.0525 .  
## PH_D         1.748e-02  7.462e-03   2.342   0.0213 *  
## GRADUAT      1.808e-02  3.432e-03   5.268 9.03e-07 ***
## OUT_STAT    -5.027e-05  2.389e-05  -2.104   0.0381 *  
## BOOK         1.044e-03  4.510e-04   2.315   0.0228 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4225 on 92 degrees of freedom
##   (25 пропущенных наблюдений удалены)
## Multiple R-squared:  0.4772, Adjusted R-squared:  0.4488 
## F-statistic:  16.8 on 5 and 92 DF,  p-value: 9.211e-12
AIC(model_pub_nono)
## [1] 117.0556

Комментарий: эта модель в сравнении с model_pub_nono имеет более высокий p-value и более низкий R^2.

plot(model_pub_full, 1:5)

4.8.7 Поиск outliers по Residuals vs. Deleted Residuals

plot(residuals(model_pub_full), rstudent(model_pub_full))

Комментарий: особых выбросов нет.

influence.measures(model_pub_full)
## Influence measures of
##   lm(formula = log_NEW10 ~ GRADUAT + PH_D + BOOK + OUT_STAT, data = pub_data) :
## 
##        dfb.1_  dfb.GRAD  dfb.PH_D  dfb.BOOK  dfb.OUT_    dffit cov.r   cook.d
## 1   -0.042131 -0.133457  3.02e-02  0.016099  1.27e-01 -0.17853 1.047 6.39e-03
## 3    0.103254  0.106100 -1.30e-01 -0.039284  4.56e-02 -0.20747 1.147 8.67e-03
## 5    0.025264  0.031849 -1.62e-02 -0.037843 -1.18e-02 -0.07080 1.072 1.01e-03
## 6    0.042930 -0.026515 -3.75e-02  0.020769  2.19e-03  0.07699 1.074 1.20e-03
## 7   -0.001421 -0.001482  2.05e-03 -0.000293  3.65e-05  0.00295 1.083 1.76e-06
## 8    0.256242 -0.005185 -1.43e-01 -0.095065 -1.30e-01  0.30174 1.034 1.81e-02
## 9   -0.223857  0.071406  9.59e-02  0.076188  2.17e-01  0.40619 0.932 3.22e-02
## 11  -0.350149 -0.136018  1.30e-01  0.312988  3.54e-01  0.52495 1.037 5.43e-02
## 12  -0.121942  0.075940  3.86e-02  0.138385 -2.94e-02  0.21075 1.137 8.94e-03
## 13  -0.215844 -0.020160  1.53e-01  0.134773 -1.49e-03  0.27751 1.055 1.54e-02
## 15  -0.321287 -0.077636  2.49e-01  0.025996  2.44e-01  0.44045 0.919 3.77e-02
## 16  -0.250015 -0.049278  2.92e-01 -0.032726 -4.84e-02  0.38214 0.833 2.80e-02
## 17  -0.081306  0.044052  6.07e-02  0.007316 -2.62e-02 -0.11473 1.077 2.65e-03
## 18   0.003194 -0.000841  1.97e-02 -0.041553 -4.86e-03  0.05185 1.085 5.43e-04
## 19   0.009679  0.007517 -7.39e-03  0.003912 -2.18e-02 -0.02766 1.122 1.55e-04
## 20   0.141867  0.108893 -1.32e-01  0.020169 -1.07e-01  0.23226 0.956 1.06e-02
## 21  -0.012381  0.071752  1.14e-02  0.006585 -6.69e-02  0.12529 1.036 3.15e-03
## 22   0.036851 -0.032196  3.27e-02 -0.086552 -3.76e-02  0.12703 1.053 3.24e-03
## 23  -0.112694  0.134024  1.63e-02  0.188038 -1.08e-01  0.30489 1.045 1.85e-02
## 24   0.134601  0.381840 -1.22e-01 -0.241573 -1.15e-01 -0.53441 0.843 5.46e-02
## 25   0.014454  0.059799  3.86e-02 -0.077258 -1.06e-01  0.13967 1.084 3.93e-03
## 26   0.006213  0.014783 -3.94e-03  0.011396 -3.38e-02  0.04477 1.100 4.05e-04
## 27  -0.003111  0.012714  4.84e-02 -0.119863  2.89e-03 -0.13920 1.130 3.91e-03
## 28  -0.224463 -0.046283  1.85e-01  0.043128  1.35e-02 -0.27959 0.924 1.53e-02
## 29  -0.199615 -0.064122  1.55e-01  0.074507  2.84e-02 -0.22725 1.057 1.03e-02
## 31   0.024836  0.118218  9.65e-05 -0.070988 -8.35e-02  0.15237 1.108 4.68e-03
## 32   0.010801 -0.068147 -4.39e-04  0.018210  2.63e-02  0.08799 1.083 1.56e-03
## 33   0.546814  0.131824 -4.82e-01 -0.062290 -6.47e-02  0.58654 0.888 6.62e-02
## 34  -0.017402 -0.014653  2.72e-02 -0.010574 -9.31e-03 -0.03488 1.106 2.46e-04
## 35   0.000434 -0.007105  1.32e-03  0.002973 -6.78e-03 -0.02361 1.072 1.13e-04
## 36  -0.030540 -0.044774  4.61e-02 -0.026909  1.86e-02 -0.07377 1.075 1.10e-03
## 38   0.038315  0.014707 -3.98e-02  0.010623 -9.85e-03  0.05040 1.086 5.13e-04
## 39   0.031172 -0.009854 -1.04e-02 -0.062093  2.45e-02 -0.11232 1.044 2.53e-03
## 41  -0.001574  0.010763 -4.89e-03  0.005513 -4.69e-04 -0.01433 1.110 4.15e-05
## 42   0.033069 -0.087660  1.88e-02 -0.011504 -2.87e-02  0.15290 1.040 4.69e-03
## 44   0.043104  0.028398 -9.65e-02  0.086072  1.48e-02 -0.12728 1.088 3.26e-03
## 45   0.027188 -0.009637 -5.15e-02  0.053262 -1.60e-03 -0.11268 1.036 2.55e-03
## 46   0.032962 -0.006529 -6.83e-03 -0.044942  8.94e-03  0.07515 1.064 1.14e-03
## 47  -0.130253  0.121583  6.48e-02  0.136264 -2.07e-01 -0.38168 0.845 2.80e-02
## 48  -0.000454 -0.146705 -1.18e-01  0.270710  8.78e-02 -0.40760 0.809 3.17e-02
## 49   0.091508 -0.006039 -1.23e-01  0.115218 -6.34e-02 -0.22342 1.046 9.98e-03
## 51   0.020780  0.007590 -8.55e-03 -0.017144 -1.30e-02  0.03116 1.082 1.96e-04
## 52   0.005839 -0.001717 -4.46e-03 -0.003079  4.69e-03  0.01013 1.100 2.07e-05
## 53  -0.002512 -0.001853  1.40e-03  0.000344  3.86e-03 -0.00499 1.094 5.04e-06
## 54   0.017187  0.003092 -5.29e-03 -0.012185 -1.72e-02  0.02580 1.111 1.35e-04
## 55  -0.007263 -0.006925  6.30e-03  0.000973  1.36e-02  0.03071 1.066 1.91e-04
## 56  -0.102366 -0.142662  4.14e-02  0.142251  1.49e-01  0.25798 0.993 1.32e-02
## 58   0.071336  0.007929 -5.81e-02 -0.012835 -6.21e-03  0.07724 1.110 1.20e-03
## 59  -0.020325  0.005678  1.20e-02  0.003531  7.44e-03 -0.02635 1.101 1.40e-04
## 61   0.026528  0.015272  3.82e-02 -0.084232 -8.32e-02  0.12538 1.081 3.17e-03
## 63   0.087394 -0.059968  1.25e-02 -0.090035 -1.38e-01 -0.24911 1.013 1.23e-02
## 64  -0.170093 -0.191554  2.91e-01 -0.103360 -5.94e-02  0.34442 1.146 2.38e-02
## 66   0.247524 -0.418324 -1.60e-01 -0.152779  3.84e-01 -0.69659 0.693 8.92e-02
## 67   0.034292  0.137055 -7.46e-02  0.059278 -8.08e-02  0.16454 1.130 5.46e-03
## 68   0.026783 -0.003626 -5.07e-02  0.024273  4.16e-02 -0.08292 1.068 1.39e-03
## 69  -0.019209 -0.001041  2.31e-02 -0.007225 -4.35e-04  0.03438 1.075 2.39e-04
## 70   0.012013  0.239901 -9.68e-03 -0.078056 -1.17e-01  0.28916 1.023 1.66e-02
## 71  -0.092085 -0.005378  1.14e-01 -0.044926 -4.20e-02 -0.14123 1.070 4.01e-03
## 72   0.093203  0.066536 -2.64e-01  0.242285  1.34e-01 -0.34264 1.152 2.35e-02
## 73  -0.158831 -0.232355  1.86e-01 -0.033745  1.25e-01 -0.32859 0.894 2.10e-02
## 74  -0.110318  0.092402  9.43e-03  0.127219 -2.19e-02 -0.23596 0.963 1.10e-02
## 75  -0.062227 -0.055431  5.68e-02  0.012618  2.93e-02 -0.09321 1.070 1.75e-03
## 77   0.002545  0.003570 -4.48e-03  0.002996 -2.57e-03 -0.00797 1.079 1.28e-05
## 79   0.000397 -0.000864  1.38e-02 -0.043803  1.39e-02 -0.05684 1.119 6.53e-04
## 80  -0.003717 -0.005744 -5.30e-04  0.022151 -1.00e-02  0.03494 1.106 2.47e-04
## 83   0.030763 -0.027421 -4.54e-02  0.016224  7.42e-02  0.10076 1.086 2.05e-03
## 84   0.051012  0.143654 -5.12e-02 -0.013781 -1.43e-01 -0.19210 1.059 7.40e-03
## 85   0.023201  0.012497 -5.37e-02  0.070858 -2.02e-02 -0.09385 1.147 1.78e-03
## 86  -0.045718  0.003991  1.31e-02  0.126471 -1.41e-01 -0.29579 0.921 1.71e-02
## 87  -0.009540 -0.062359  6.97e-02 -0.098075  2.55e-03 -0.12517 1.147 3.16e-03
## 88   0.014200  0.014410 -2.16e-04 -0.027171 -1.11e-02  0.04109 1.078 3.41e-04
## 90   0.008670 -0.012445  1.30e-02 -0.071909  3.83e-02 -0.10831 1.080 2.36e-03
## 91   0.029025 -0.048617  7.56e-02 -0.120907 -9.13e-02  0.17578 1.144 6.23e-03
## 92   0.430023 -0.528612 -4.13e-01  0.206695  3.71e-01  0.83959 1.084 1.38e-01
## 93   0.045189 -0.017476  1.07e-02 -0.050952 -6.48e-02  0.10855 1.083 2.37e-03
## 96   0.196067 -0.069055 -8.44e-02 -0.082115 -9.60e-02  0.26129 1.081 1.37e-02
## 97  -0.009366  0.027867  1.62e-02 -0.002606 -4.52e-02  0.05736 1.107 6.65e-04
## 98   0.027003  0.036114 -7.06e-03 -0.057481 -1.77e-02 -0.07484 1.166 1.13e-03
## 100  0.075991  0.026069  1.59e-02  0.035124 -3.98e-01 -0.48472 1.107 4.67e-02
## 101 -0.028942  0.123097 -2.87e-02  0.017744  3.96e-02  0.18388 1.111 6.80e-03
## 102 -0.004197  0.007296  5.48e-03 -0.003717 -9.50e-03 -0.01427 1.096 4.12e-05
## 103  0.007016 -0.020480 -1.38e-03 -0.009095  2.28e-02  0.03411 1.105 2.35e-04
## 104  0.123240  0.179060 -1.41e-01  0.026554 -1.63e-01 -0.24651 1.098 1.22e-02
## 105  0.023884 -0.063654  8.23e-02 -0.146901 -7.63e-02 -0.20001 1.083 8.03e-03
## 106  0.055488 -0.001506 -4.34e-02 -0.022109  5.94e-04 -0.07177 1.090 1.04e-03
## 107  0.004810  0.000641  1.61e-02 -0.041155 -1.16e-02 -0.04542 1.121 4.17e-04
## 109  0.015535 -0.011063 -2.23e-02  0.019857  3.92e-03 -0.03938 1.093 3.13e-04
## 110 -0.045418  0.099942  3.16e-02 -0.007424 -8.10e-02 -0.14449 1.071 4.20e-03
## 111  0.003248  0.005327 -1.45e-03 -0.004667 -4.85e-03 -0.00846 1.105 1.45e-05
## 113  0.037474 -0.129469  3.64e-03 -0.054592  7.57e-02 -0.17693 1.065 6.28e-03
## 114  0.223708  0.055697 -2.56e-01 -0.035107  2.07e-01  0.42656 0.901 3.53e-02
## 115 -0.039680  0.195298  2.60e-03 -0.087338  5.77e-02  0.32221 1.088 2.07e-02
## 116 -0.012937 -0.026990  1.32e-02  0.008373  9.61e-03 -0.03206 1.140 2.08e-04
## 117 -0.137025  0.039881  4.03e-02 -0.106701  4.40e-01  0.61275 1.075 7.40e-02
## 118  0.001980 -0.000154 -1.45e-03 -0.002237  1.56e-03 -0.00497 1.084 4.99e-06
## 119 -0.195111  0.041427  1.92e-01 -0.064464 -1.77e-02 -0.25500 1.043 1.30e-02
## 122 -0.049622 -0.115513  9.73e-02 -0.034789  2.04e-02  0.14702 1.104 4.35e-03
## 123  0.002463  0.058149  1.39e-02 -0.020111 -7.73e-02  0.09267 1.112 1.73e-03
##        hat inf
## 1   0.0367    
## 3   0.0988    
## 5   0.0252    
## 6   0.0278    
## 7   0.0254    
## 8   0.0578    
## 9   0.0482    
## 11  0.1039    
## 12  0.0931    
## 13  0.0613    
## 15  0.0513    
## 16  0.0283   *
## 17  0.0378    
## 18  0.0312    
## 19  0.0595    
## 20  0.0232    
## 21  0.0210    
## 22  0.0278    
## 23  0.0626    
## 24  0.0515    
## 25  0.0471    
## 26  0.0428    
## 27  0.0772    
## 28  0.0256    
## 29  0.0513    
## 31  0.0647    
## 32  0.0365    
## 33  0.0693    
## 34  0.0470    
## 35  0.0171    
## 36  0.0278    
## 38  0.0325    
## 39  0.0210    
## 41  0.0491    
## 42  0.0283    
## 44  0.0472    
## 45  0.0181    
## 46  0.0214    
## 47  0.0296    
## 48  0.0288   *
## 49  0.0460    
## 51  0.0267    
## 52  0.0403    
## 53  0.0354    
## 54  0.0508    
## 55  0.0137    
## 56  0.0355    
## 58  0.0543    
## 59  0.0418    
## 61  0.0420    
## 63  0.0393    
## 64  0.1210    
## 66  0.0505   *
## 67  0.0814    
## 68  0.0257    
## 69  0.0209    
## 70  0.0510    
## 71  0.0391    
## 72  0.1239    
## 73  0.0286    
## 74  0.0248    
## 75  0.0290    
## 77  0.0218    
## 79  0.0593    
## 80  0.0472    
## 83  0.0403    
## 84  0.0447    
## 85  0.0845    
## 86  0.0277    
## 87  0.0877    
## 88  0.0242    
## 90  0.0379    
## 91  0.0922    
## 92  0.1798   *
## 93  0.0399    
## 96  0.0703    
## 97  0.0495    
## 98  0.0974   *
## 100 0.1264    
## 101 0.0723    
## 102 0.0374    
## 103 0.0464    
## 104 0.0760    
## 105 0.0584    
## 106 0.0386    
## 107 0.0597    
## 109 0.0362    
## 110 0.0403    
## 111 0.0452    
## 113 0.0442    
## 114 0.0452    
## 115 0.0860    
## 116 0.0745    
## 117 0.1361    
## 118 0.0261    
## 119 0.0515    
## 122 0.0610    
## 123 0.0579

4.8.8 Поиск outliers по Cook и Mahalanobis

4.8.8.1 Cook

# Рассчитаем Cook's distance для модели с логарифмами
cooks_vals <- cooks.distance(model_pub_full)
cook_df <- data.frame(
  obs  = row.names(model_pub_full$model),
  cook = cooks_vals
) %>%
  arrange(desc(cook))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(cook_df))


plot(
  x, cook_df$cook,
  type  = "h",
  main  = "Cook's Distance",
  xlab  = "Номер наблюдения (из исходного df)",
  ylab  = "Cook's distance",
  xaxt  = "n"                # убираем автоматическую ось X
)
axis(side = 1, at = x, labels = cook_df$obs, las = 2, cex.axis = 0.8)
abline(h = 1, col = "red", lty = 2)

Комментарий: высокие значения у 29 и 38.

4.8.8.2 Mahalanobis

# 7.2. Расстояние Махаланобиса (Mahalanobis Distance)
# Расстояние Махаланобиса в пространстве предикторов (регрессоров)
mahalanobis_distance <- mahalanobis(model.matrix(model_pub_full)[,-1],
                                     center = colMeans(model.matrix(model_pub_full)[,-1]),
                                     cov = cov(model.matrix(model_pub_full)[,-1]))
mah_df <- data.frame(
  obs  = row.names(model_pub_full$model),
  mah = mahalanobis_distance
) %>%
  arrange(desc(mah))
# Создадим вектор x для индексов (уже упорядоченных)
x <- seq_len(nrow(mah_df))
alpha<-0.05

# Scatter plot для расстояния Махаланобиса
plot(
  x, mah_df$mah,
  type  = "h",
  main  = "Mahalanobis's Distance",
  xlab  = "Номер наблюдения (из исходного df)",
  ylab  = "Mahalanobis's distance",
  xaxt  = "n"                
)
axis(side = 1, at = x, labels = mah_df$obs, las = 2, cex.axis = 0.8)
abline(h = qchisq(p = 1 - alpha, df = nrow(mah_df)), col = "red", lty = 2)

4.8.9 Итог: прогнозирование

set.seed(456) # Другое зерно для генерации новых значений
max_GRADUAT<-max(pub_data$GRADUAT, na.rm=TRUE)
max_BOOK<-max(pub_data$BOOK, na.rm=TRUE)
max_OUT_STAT<-max(pub_data$OUT_STAT, na.rm=TRUE)
max_PH_D<-max(pub_data$PH_D, na.rm=TRUE)
new_data <- data.frame(
  GRADUAT = max_GRADUAT - 1:8,
  BOOK = max_BOOK - 1:8,
  OUT_STAT = max_OUT_STAT - seq(100, 800, by=100),
  PH_D = max_PH_D - 1:8
)
# Прогнозируем значение NEW10
predicted_values <- predict(model_pub_full, newdata = new_data)
cat("Прогнозируемые значения NEW10:\n")
## Прогнозируемые значения NEW10:
print(predicted_values)
##        1        2        3        4        5        6        7        8 
## 4.330459 4.291791 4.253122 4.214454 4.175785 4.137117 4.098448 4.059780
# 6. Доверительные и предсказательные интервалы

# Доверительный интервал (confidence interval) для среднего значения NEW10
confidence_intervals <- predict(model_pub_full, newdata = new_data, interval = "confidence", level = 0.95)
print("Доверительные интервалы (95%):\n")
## [1] "Доверительные интервалы (95%):\n"
print(confidence_intervals)
##        fit      lwr      upr
## 1 4.330459 3.917940 4.742979
## 2 4.291791 3.883437 4.700145
## 3 4.253122 3.848577 4.657668
## 4 4.214454 3.813349 4.615558
## 5 4.175785 3.777744 4.573826
## 6 4.137117 3.741753 4.532480
## 7 4.098448 3.705369 4.491527
## 8 4.059780 3.668584 4.450975
# Предсказательный интервал (prediction interval) для отдельного наблюдения NEW10
prediction_intervals <- predict(model_pub_full, newdata = new_data, interval = "prediction", level = 0.95)
print("Предсказательные интервалы (95%):\n")
## [1] "Предсказательные интервалы (95%):\n"
print(prediction_intervals)
##        fit      lwr      upr
## 1 4.330459 3.384023 5.276896
## 2 4.291791 3.347162 5.236419
## 3 4.253122 3.310134 5.196111
## 4 4.214454 3.272936 5.155971
## 5 4.175785 3.235569 5.116001
## 6 4.137117 3.198031 5.076202
## 7 4.098448 3.160322 5.036574
## 8 4.059780 3.122441 4.997118
# Создаем data.frame для прогнозов и интервалов
predictions_df <- data.frame(
  GRADUAT = new_data$GRADUAT,
  OUT_STAT = new_data$OUT_STAT,
  BOOK = new_data$BOOK,
  PH_D = new_data$PH_D,
  Predicted = predicted_values,
  Lower_Conf = confidence_intervals[, "lwr"],
  Upper_Conf = confidence_intervals[, "upr"],
  Lower_Pred = prediction_intervals[, "lwr"],
  Upper_Pred = prediction_intervals[, "upr"]
)
predictions_df
# Визуализация (только для демонстрации, нужен более сложный график для нескольких предикторов)
ggplot(predictions_df, aes(x = 1:nrow(predictions_df), y = Predicted)) + # Условная ось x, т.к. несколько предикторов
  geom_point() +
  geom_errorbar(aes(ymin = Lower_Conf, ymax = Upper_Conf), color = "blue", width = 0.2) + # Доверительные интервалы
  geom_errorbar(aes(ymin = Lower_Pred, ymax = Upper_Pred), color = "red", width = 0.2) + # Предсказательные интервалы
  labs(title = "Прогнозы NEW10 с доверительными и предсказательными интервалами",
       x = "Прогноз (номер)",
       y = "NEW10") +
  scale_x_continuous(breaks = 1:nrow(predictions_df)) + # Подписи для оси x
  theme_bw()

  # annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.9, label = "Синие: Доверительные интервалы", color = "blue", hjust = 0) +
  # annotate("text", x = 1, y = max(predictions_df$Upper_Pred) * 0.8, label = "Красные: Предсказательные интервалы", color = "red", hjust = 0)